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Coherent  Lidar  System  for  High  Resolution  Measurement  of 
Atmospheric  Wind  Velocity  and  Water  Vapor  Fields 


I.  SUMMARY 
1.1  OVERVIEW 

The  goal  of  this  study  is  to  evaluate  the  capability  for  a  scanning  coherent  pulsed  laser 
radar/lidar,  to  achieve  a  high  resolution  mapping  of  the  three  dimensional  structure  of 
wind  velocity  and  water  vapor  in  the  turbulent  planetary  boundary  layer.  This 
information  is  needed  to  validate  models  of  turbulent  dispersion  in  the  atmospheric 
boundary  layer  under  a  variety  of  stability  and  terrain  conditions.  It  is  generally 
recognized  that  intermittent  or  sporadic  coherent  events  in  the  boundary  layer  may  play 
a  vital  part  in  determining  the  overall  transport  mechanisms,  both  of  species  as  well  as 
momentum  and  energy.  Large  scale  3D  computations  of  boundary  layer  dynamics  are 
playing  an  increasing  role  in  the  development  of  the  understanding  of  these  transport 
processes.  It  is  of  particular  interest  to  obtain  quantitative  observations  of  the  eddy 
transport  mechanics  at  the  grid  scales  of  these  computations  (several  tens  of  meters) 
so  that  comparisons  can  be  made  and  models  of  the  transport  can  be  validated. 

Since  its  inception  in  1984  Coherent  Technologies,  Inc.  (CTI)  has  been  developing 
scanning  solid-state  coherent  lidars  that  have  the  ability  to  probe  the  atmospheric 
boundary  layer  at  scales  relevant  to  the  description  of  atmospheric  boundary  layer  wind 
fields.  The  basic  objective  of  the  present  study  is  to  determine  the  feasibility  and  to 
construct  a  conceptual  design  of  a  lidar  system  for  measuring  both  the  turbulence  itself 
(wind  velocity  fluctuations)  and  its  consequence  (transport)  to  meet  the  above 
objectives.  For  this  purpose  the  study  has  been  divided  into  four  tasks: 

1.  Generate  realistic  examples  of  wind  velocity  and  water  vapor 
turbulent  or  fluctuation  fields  and  simulate  the  properties  of  a 
scanning  lidar  diagnostic  system  to  make  Doppler  measurement  of 
line-of-sight  velocity  and  DIAL  (Differential  Absorption  Lidar) 
measurements  of  water  vapor  concentration. 

2.  Predict  performance  and  measurement  precisions  to  be  expected 
for  line  of  sight  velocity  and  concentrations  as  a  function  of  space 
and  time  resolution,  scan  coverage,  atmospheric  environment,  and 
lidar  parameters. 

3.  Demonstrate  velocity  measurement  capability  in  the  boundary  layer 
and  validate  basic  parameters  of  the  prediction  models  and 
measurement  concepts  using  available  Doppler  lidar  scanning 
systems. 

4.  Construct  a  conceptual  design  of  a  prototype  system  for  turbulence 
measurements  in  the  planetary  boundary  layer. 


2 


There  are  basically  three  lidar  measurement  and  data  analysis  concepts  that  can  be 
considered  for  measuring  properties  of  the  turbulent  atmospheric  boundary  layer  In 
order  of  increasing  complexity  these  are: 

1.  Conventional  sampling  at  one  or  more  space  points  from  one  or 
more  view  directions  of  the  time  history  of  the  local  radial  (line-of- 
sight)  wind  followed  by  a  multidimensional  spectral  or  conditional 
sampling  analysis  to  construct  the  statistical  properties  of  the  local 
wind  field. 

2.  Measurement  of  the  spatial  distribution  of  the  radial  wind  component 
over  a  3  dimensional  measurement  field  followed  by  a  data  analysis 
that  imposes  continuity,  persistence,  spatial  uniformity  or  other 
constraints  to  restrict  the  form  of  allowed  solutions  and  to  fit  these 
forms  to  the  data. 

3.  Incorporation  of  the  radial  wind  components  measured  by  the  lidar 
as  a  function  of  space  and  time  into  a  full-blown  3D  hydrodynamic 
calculation  of  the  boundary  layer  flow. 

A  major  result  of  the  present  study  is  the  conclusion  that  lidar  and  signal  processing 
capabilities  that  have  recently  become  available  allow  realistic  implementation  of  the 
third  concept.  The  first  two  concepts  have  been  extensively  used  in  the  past,  require 
limited  sampling  of  the  flow  environment,  and  can  be  implemented  with  limited 
scanning  and  signal  processing  tools.  The  third  requires  a  full  3D  scanning  capability 
together  with  a  signal  processing  that  can  retrieve  the  3D  distribution  of  the  vector  wind 
field.  The  second  method,  which  includes  such  concepts  as  VAD  and  cloud  tracking, 
is  essentially  a  derivative  of  the  third  with  idealized  assumptions  introduced  to  limit  the 
measurement  and  processing  requirements  (such  as  horizontally  homogeneous  winds 
and/or  ’frozen’  turbulence).  The  first  method  is  a  statistical  treatment  of  the  wind  field 
measurements  designed  to  extract  stationary  measures  of  the  flow  and  can  be  applied 
directly  to  the  measurements  in  real  time  or  to  recorded  data.  Given  that  the  third 
method  is  possible  it  provides  data  that  can  be  processed  after  the  measurement  using 
such  statistical  techniques. 

This  report  will  emphasize  the  third  approach.  To  retrieve  both  vector  winds  and  water 
vapor  concentration  at  small  turbulence  scales,  we  have  found  it  necessary  to  develop 
two  new  signal  processing  concepts.  First,  a  processing  concept  has  been  developed 
to  retheve  three  component  vector  winds  from  the  one  component  line-of-sight  velocity 
measurements  by  the  imposition  of  known  hydrodynamic  constraints  (mass 
conservation  and  momentum  conservation).  This  procedure  is  basically  an  extension 
of  a  well-known  method  used  to  retrieve  vector  winds  from  two  station  radar  data  (Dual 
Doppler  method).1  For  application  to  lidar  we  have  extended  this  method  to  apply  to 
single  station  data  by  structuring  the  signal  processing  analysis  in  a  Bayesian 
estimation  framework  that  merges  the  measurements  with  a  complete  description  of  the 
equations  governing  the  boundary  layer  dynamics. 


*R.J.  Doviak  and  D.S.  Zrnic,  "Doppler  Radar  and  Weather  Observations,”  Chapt.  9,  Academic  Press 
(1984) 
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A  second  measurement  and  processing  concept  has  been  developed  to  allow  Doppler 
velocity  data  and  DIAL  measurements  of  water  vapor  concentration  to  be  combined 
and  interpreted  in  terms  of  small  scale  correlations  and  transport.  The  major  problem 
with  DIAL  for  turbulence  studies  is  its  limited  capability  for  sensing  small  scale 
structures.  DIAL  is  an  indirect  method  for  sensing  species  concentrations  that  relies  on 
spatial  and  spectral  differencing  to  interpret  absorption  features  of  the  reflected  light  in 
terms  of  local  species  concentrations.  This  differencing  at  small  spatial  scales 
represents  a  major  source  of  noise  at  the  weak  reflection  levels  characteristic  of 
atmosphere  aerosol  backscatter  and  effectively  precludes  the  making  of  accurate 
measurements  at  turbulence  scales.  The  technique  proposed  here  utililizes  DIAL  to 
establish  a  correlation  at  large  scales  between  water  vapor  density  fluctuations  and 
fluctuations  of  the  backscattering  coefficient  due  to  aerosols.  Both  are  passively 
advected  by  the  turbulent  wind  field  and  similar  correlation  is  to  be  expected  at  both 
large  and  small  scales.  Since  the  fluctuations  of  aerosol  density  are  much  more  easily 
measured  at  small  scales  than  are  the  absorption  fluctuations,  an  effective  method  of 
extending  DIAL  to  small  scales  is  to  measure  both  absorption  (using  DIAL)  and 
backscatter  at  large  scales,  establish  a  correlation  coefficient,  and  use  the  correlation 
coefficient  to  infer  small  scale  water  vapor  fluctuations  from  similar  fluctuations  in  the 
measured  backscatter. 

Evaluation  of  the  expected  performance  and  simulation  of  the  data  processing  have 
been  carried  out  to  validate  the  concepts  and  determine  the  system  measurement 
capability.  The  results  of  the  analyses  are  described  in  the  following  section. 


1.2  RESULTS  SUMMARY 

1.2.1  VELOCITY  ERROR  EXPECTATIONS 

At  high  SNR  the  precision  with  which  the  velocity  can  be  estimated  in  a  single  pulse  is 
determined  by  the  width  of  the  Doppler  spectrum.  For  short  pulses  the  width  is 
determined  by  the  pulse  duration.  For  long  pulses  the  spectrum  represents  the  spread 
of  aerosol  velocities  within  the  pulse.  Figure  1.1  shows  the  dependence  of  this  width 
on  the  subgrid  turbulence  level  and  lidar  wavelength.  For  a  pulsed  system  in  the 
absence  of  turbulence,  the  product  of  expected  velocity  error  standard  deviation  (per 
pulse)  and  pulse  length  is  a  constant. 

For  atmospheric  turbulence  measurements  we  have  assumed  that  a  velocity  precision 
of  the  order  of  1.0  meter/second  or  better  per  pulse  is  adequate  to  monitor  the 
turbulence  when  many  pulses  sample  the  3D  observational  volume.  Greater  velocity 
precision  can  be  achieved  by  averaging  pulses  at  the  expense  of  decreased  spatial 
resolution  or  coverage.  Figure  1.1  shows  that  this  requirement  limits  the  achievable 
spatial  range  resolution  to  about  50  meters  at  2  microns  wavelength.  For  a  given 
velocity  precision,  the  spatial  resolution  is  proportional  to  wavelength  in  the  absence  of 
sub  grid  turbulence  (i.e.,  250  meters  at  10  microns  and  25  meters  at  1  micron).  In 
general,  shorter  wavelengths  are  preferred  for  a  pulsed  system  when  both  high  velocity 
and  spatial  resolution  are  required.  In  numerical  simulations  of  the  velocity  retrieval  we 
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Figure  1.1  Spectrum  width.  The  Doppler  spectrum  width  can  be  minimized  in  the 
presence  of  turbulence  by  appropiate  choice  of  pulse  length.  Smaller 
widths  can  be  achieved  at  smaller  wavelengths 


have  found  that  maintaining  a  high  lateral  resolution  also  assists  in  maintaining  a  high 
range  resolution  as  a  result  of  the  coupling  to  the  fluid  dynamic  constraints  of  the  flow. 

In  the  presence  of  turbulence,  the  expected  velocity  error  shows  a  minimum  as  a 
function  of  pulse  length  at  the  point  where  the  spectral  broadening  due  to  the  finite 
pulse  length  is  comparable  to  that  created  by  the  spread  of  aerosol  radial  velocities 
within  the  pulse.  For  light  to  moderate  turbulence  this  occurs  at  a  pulse  length  of  order 
100  meters  at  2  microns.  The  achievable  velocity  resolution  is  then  relatively 
insensitive  to  pulse  lengths  over  the  range  50  to  150  meters. 

In  general  we  expect  that  achievable  spatial  resolution  in  the  direction  of  look  for  a 
pulsed  coherent  Doppler  lidar  designed  to  probe  atmospheric  turbulence  will  be  limited 
to  about  25  to  50  meters  at  1  and  2  microns  respectively.  Some  benefit  may  result 
from  high  density  angular  sampling  and/or  shorter  pulses  with  multiple  pulse  averaging. 
On  the  basis  of  these  considerations  we  have  assumed  a  nominal  design  requirement 
of  25  to  50  meters  lateral  separation  of  the  lines  of  sight.  This  will  lead  to  the 
specification  of  at  least  300  Hz  prf  with  at  least  1  watt  mean  power  to  probe  a  volume 
1  km  x  1  km  in  lateral  extent  and  1  to  3  km  long  with  a  3  second  update  period.  Both  a 
one  micron  and  a  two  micron  system  are  candidates.  At  ten  microns  focusing  with  long 
or  CW  pulses  is  needed  to  achieve  adequate  range  resolution.  This  will  severely  limit 
the  maximum  allowable  range.  In  general  one  micron  has  best  resolution  under  ideal 
conditions  but  is  not  eyesafe  and  is  more  readily  degraded  by  refractive  turbulence. 


1.2.2.  SIMULATIONS 

Several  simulation  models  are  available  at  CTI  to  characterize  the  lidar  reponse  to  both 
localized  wind  and  aerosol  patterns.  These  models  provide  estimates  of  the 
measurement  precision  expected  as  a  function  of  lidar  system  design  and  atmospheric 
state.  Simulations  of  candidate  turbulence  wind  fields  have  been  generated  in  two 
dimensions  and  have  been  used  in  a  simulation  of  the  measurement  concept.  In  these 
simulations  a  two  dimensional  Kalman  filter  accepts  radial  velocity  measures  made  at 
arbitrary  space  and  time  coordinates  in  a  vertical  plane  and  incorporates  these  into  a 
numerical  2D  incompressible  hydro  model  of  the  boundary  layer.  The  2D  code  is  not 
appropriate  for  the  real  atmospheric  boundary  layer  but  is  being  used  to  test  the  data 
assimilation  procedures.  A  full  3D  code  will  be  required  for  the  actual  assimilation  and 
nowcasting  of  real  data.  In  the  2D  analysis  each  line  of  sight  is  intended  to  represent  a 
plane  of  3D  data. 

Figures  1.2  and  1.3  show  examples  of  the  operation  and  type  of  output  of  this  analysis. 
Measurements  of  the  radial  wind  are  assumed  to  be  continuously  collected  along 
several  lines  of  sight.  In  the  simulation  these  wind  fields  were  generated  by  a  2D 
hydrodynamic  boundary  layer  code.  These  simulated  measurements  are  fed  to  a 
Kalman  filter  which  continuously  updates  a  similar  2D  hydro  code  describing  the  vector 
wind  field  in  a  vertical  plane  containing  the  measurement  points.  The  difference 
between  the  measured  data  and  the  code  predicted  data  at  the  measurement  points 
forms  a  source  of  distributed  vorticity  throughout  the  computational  volume  that  allows 
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(a) 


Figure  1.2.  Kalman  filter  simulation  to  adapt  a  2D  boundary  layer  hydro  code  to 
simulated  lidar  measurements  of  the  near  surface  line-of-sight  winds. 

(a)  Velocity  patterns  and  line-of-sight  velocity  components  at  51  seconds 
The  lower  panels  compare  the  simulated  measurements(solid  line)  with 
the  current  estimate  (dots)  provided  by  the  nowcast  output  of  the  hydro 
code.  Panels  marked  with  an  asterisk  identify  participating  lines  of  sight. 

(b)  Evolution  of  the  estimated  2D  field  at  101  seconds.  At  this  time  the 
predicted  data  are  in  essential  agreement  with  the  measured  data. 
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(C) 


Figure  1.3.  Kalman  filter  simulation  to  adapt  a  2D  boundary  layer  hydro  code  to 
simulated  lidar  measurements  of  line-of-sight  winds. 

(a)  True  velocity  patterns  at  time  101  seconds 

(b)  Velocity  patterns  and  line-of-sight  velocity  components  inferred  from 
three  lines  of  sight. 

(c)  Velocity  patterns  and  line-of-sight  velocity  components  inferred  from 
six  lines  of  sight. 
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the  atmosphere  to  adapt  to  the  measurements.  The  algorithms  and  implementation  of 
this  procedure  is  described  in  detail  in  the  Section  3.1. 

Figure  1.2  shows  the  adaption  to  a  wind  field  created  by  a  simple  vortex  pair 
representation  of  a  wind  gust.  Simulated  data  has  been  acquired  along  three  lines  of 
sight.  In  the  reconstruction,  the  inferred  atmosphere  is  assumed  to  be  initially 
motionless  at  time  0.  The  Kalman  filter  then  generates  and  advects  vorticity  to  allow 
the  atmosphere  to  adapt  to  the  measurements.  In  this  example  the  vector  velocity  field 
becomes  well  reconstructed  throughout  the  region  displayed  within  100  seconds  after 
initation  of  the  adaptation. 

Reconstructions  of  a  more  complex  wind  field  are  shown  in  Figure  1.3.  The  true 
velocity  field  at  time  100  seconds  is  shown  in  panel  a  and  reconstructions  using  3  and 
6  lines  of  sight  in  panels  b  and  c.  In  this  case  the  reconstruction  with  3  lines  of  sight  is 
less  satisfactory  than  in  the  prior  example  even  though  the  deduced  field  at  the 
measurement  points  is  very  accurate.  With  6  lines  of  sight  the  inferred  field  is  well 
predicted  throughout  the  measurement  field. 

As  is  apparent  from  these  examples  the  inferred  field  is  one  possible  solution,  but  not 
the  only  one.  When  only  a  small  number  of  lines  of  sight  are  measured,  there  are 
potentially  many  possible  atmospheric  flows  that  can  reproduce  the  measurements.  As 
demonstrated  in  Figure  1.3  fine  gridding  of  the  lines  of  sight  can  reduce  the  level  of 
ambiguity  and  improve  the  quality  of  the  inferred  field. 

These  calculations  demonstrate  that  a  hydrodynamically  constrained  analysis  can  be 
used  to  generate  vector  winds  from  measurements  of  only  one  component. 


1.3  WATER  VAPOR  MEASUREMENT  USING  DIAL 

DIAL  meaurements  of  the  type  discussed  in  Section  1.1  require  a  high  prf  source 
(typically  300  Hz)  that  can  be  switched  between  an  absorbing  (A)  and  non-absorbing 
(N)  wavelength.  In  selecting  suitable  wavelengths  for  DIAL  measurements,  several 
issues  must  be  considered.  First,  the  non-absorbing  one  should  have  quite  a  low 
attenuation  over  the  desired  range.  Second,  the  absorbing  one  should  have  an 
attenuation  that  is  not  too  weak,  but  also  not  too  strong.  If  it  is  weak,  the  differential 
absorption  becomes  more  difficult,  while  if  it  is  too  strong,  there  may  not  be  enough 
signal  to  make  an  accurate  estimate  of  water  vapor  from  the  furthest  range.  The  A 
wavelength  must  also  be  directly  associated  with  water  absorption,  and  should  ideally 
not  be  contaminated  by  absorption  due  to  other  species. 

Figure  1.4a  shows  the  MLS  (Mid-latitude  Summer)  attenuation  spectrum  in  the  region 
of  tuning  of  TnrYAG  lasers.  There  is  a  clear  C02  overtone  band  covering  much  of  the 
region,  as  well  as  other  peaks  associated  with  water  and  N02.  Figure  1.4b  shows  the 
water  spectrum  alone.  Several  of  the  water  lines  overlap  directly  with  C02  lines  and 
are  therefore  unsuitable  for  the  present  application.  The  line  marked  with  an  arrow  and 
the  letter  A  at  2016.974  nm  (vacuum  wavelength)  has  a  peak  extinction  coefficient  of 
3.389  (per  km  roundtrip)  and  is  virtually  free  from  contamination  by  neighboring  C02 
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lines  (these  lines  contribute  <  3%  of  the  total  attenuation.  It  is  therefore  an  excellent 
choice  for  short-range  water  vapor  measurements. 

For  the  non-absorbing  wavelength,  the  best  choice  is  to  operate  at  2017.713  nm.  The 
valley  shows  a  total  absorption  of  0.1  per  km  at  that  point  (0.4  db  per  km  roundtrip). 
The  width  of  the  valley  floor  is  approximately  0.078  nm,  corresponding  to  5.8  GHz. 

Both  of  these  wavelengths  are  easily  accessible  with  an  injection-seeded  Tm:YAG 
laser.  Continuous  tuning  from  2006  to  2023  nm  has  been  demonstrated  at  CTI. 

1.4  CANDIDATE  MEASUREMENT  AND  PROCESSING  CONCEPTS 

The  ultimate  lidar  boundary  layer  sensor  would  have  the  ability  to  probe  a  three 
dimensional  volume  and  produce  a  detailed  map  of  the  line  of  sight  velocity  field  with  of 
order  20  meter  resolution  throughout  a  1  kilometer  cube  of  atmosphere  at  a  refresh  rate 
of  0.2  to  0.5  Hz.  Such  a  direct  data  sensing  system  would  operate  at  up  to  a  kilohertz 
pulse  repetition  frequency,  and  would  require  up  to  10  watts  mean  power,  and  would 
require  multiple  ground  installations  to  retrieve  the  full  vector  field. 

The  basic  premise  of  the  current  analysis  is  that  a  much  less  capable  system  in  terms 
of  hardware  can  be  expected  to  perform  the  same  basic  interrogation  task  if  combined 
with  suitable  data  processing.  The  required  data  processing  must  recognize  that  the 
'normal'  atmosphere  dynamics  are  severely  limited  by  requirements  of  continuity  and 
momentum  conservation.  Appropriate  merging  of  these  constraints  with  the  measured 
data  can  substantially  lessen  the  sensing  requirements  without  materially  restricting  the 
validity  of  the  intepretation  of  the  observation.  The  goal  of  the  present  study  is  to 
define  and  test  examples  of  such  combined  measurement  and  data  analysis  concepts 
using  available  scanning  lidar  systems. 

We  envisage  monitoring  the  radial  wind  over  multiple  vertical  or  horizontal  planes,  each 

1.5  km  wide  by  3  km  long,  centered  at  a  horizontal  range  of  3  or  4  km.  With  a  mesh 
spacing  equal  to  the  pulse  length  (50  meters)  the  advection  time  across  the  mesh  cell 
at  a  mean  wind  speed  of  10  meters/sec  is  5  seconds.  Thus,  we  will  need  an  update 
rate  of  something  like  0.2  Hz  to  follow  advecting  wind  patterns.  This  coverage  will 
require  about  30  pulses  per  plane  every  5  seconds  (a  prf  of  6  pulses  per  plane  per 
second)  and  would  allow  a  lateral  coverage  of  a  1.5  km  and  a  longtitudinal  coverage  of 
at  least  3  km.  Thus  we  envisage  monitoring  the  radial  wind  over  multiple  vertical  or 
horizontal  planes,  each  1.5  km  wide  by  3  km  long  centered  at  a  horizontal  range  of  3 
or  4  km.  With  a  mean  wind  speed  of  10  meters  per  second  this  allows  at  total 
observation  time  will  be  2  to  5  minutes  before  the  scene  would  entirely  decorrelate 
(total  of  750  to  1500  pulses  transmitted  per  plane  and  available  for  data  analysis  -  a  net 
density  of  30  range  gates  per  plane,  spread  out  in  time,  per  horizontal  resolution  cell). 

The  issue  is  how  many  measurement  planes  will  be  required.  Figures  1.5  shows  the 
expected  measurement  characteristics  for  a  candidate  measurement  that  uses 
relatively  dense  sampling  (40  planes).  Here  a  2  km  x  2  km  x  2  km  cube  centered  at  3 
km  horizontal  range  is  to  be  scanned  by  a  single  ground  lidar.  The  lidar  is  assumed  to 
operate  at  2  watts  mean  power,  with  a  320  pulse  per  sec  repetition  frequency,  and  with 
a  pulse  length  of  0.3  microseconds  (FWHM).  At  a  refresh  rate  of  0.2  Hz  this  provides  a 
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40  x  40  x  40  mesh  with  an  average  of  50  meter  resolution.  The  pulse  energy  is  6.7 
millijoules  and  is  50  meters  long. 


1.5  PRELIMINARY  VELOCITY  MEASUREMENTS  AND  PROCESSING 

It  would  be  highly  desirable  to  demonstrate  the  velocity  retrieval  process  using  existent 
lidar  equipment.  The  simulations  shown  in  Figures  1.2  to  1.3  indicate  that  a  reasonable 
reconstruction  of  the  vector  wind  field  might  be  achievable  with  as  few  as  6  scan 
planes.  This  would  correspond  to  a  36  prf,  0.25  watt  system. 


Vertical  scan  plane 

Figure  1.6.  Single  plane  measurement  geometries 

A  scanning  lidar  with  even  this  capability  is  not  presently  available  although  adequate 
systems  are  currently  under  development.  However,  two  lidars  of  lower  prf  are 
available  and  have  been  used  to  collect  data  that  can  be  used  to  demonstrate  the  data 
quality  and  some  of  the  proposed  processing  concepts. 

Two  sets  of  data  have  been  collected.  These  are  described  in  some  depth  in  Section  4 
and  are  summarized  here. 
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The  first  data  set  was  collected  using  a  flash  lamp  pumped  2  micron  laser  operating  at 
5  Hz  prf.  In  these  data  a  single,  almost  horizontal  plane  was  scanned  azimuthally  over 
about  90  degrees  and  the  data  processed  to  show  the  Doppler  signature  vs  range  at 
one  azimuth  as  well  as  the  radial  velocity  vs  range  and  azimuth.  A  second  data  set 
was  collected  using  a  flash  lamp  pumped  1  micron  laser  operating  at  10  Hz.  This  latter 
data,  taken  at  a  different  location,  scanned  a  sector  of  a  vertical  plane  in  a  nodding 
type  of  scan.  (These  pulsed  solid-state  coherent  lidar  systems  were  developed  on  U.S. 
Air  Force  contracts.  The  technical  contract  monitor  for  these  systems  is  Richard  D. 
Richmond,  Wright  Laboratories  in  Dayton,  Ohio.) 

Figures  1.7a  and  1.7b  show  two  examples  of  wind  velocity  measurements  made  with 
the  2  micron  system.  Figure  1.7a  is  a  plot  of  the  signal  frequency  content  (spectrum) 
versus  range  along  a  single  line  of  sight.  The  lidar  was  pointed  slightly  above 
horizontal  (+3°  elevation  angle)  toward  a  distant  hillside.  A  Fast  Fourier  Transform 
(FFT)  is  used  to  construct  a  frequency  spectrum  in  each  of  64  96  meter  range  gates. 
The  red  regions  indicate  the  location  in  frequency  and  range  of  the  energetic  portions 
of  the  signal  returns.  The  blue  regions  indicate  the  absence  of  signal  and  are  due  to 
the  local  oscillator  shot  noise  floor.  The  total  dynamic  range  (blue  to  red)  is  60  db.  At  a 
range  near  5  km,  an  extremely  strong  return  centered  about  0  MHz  indicates  the  distant 
hillside.  A  majority  of  the  aerosol  returns  are  centered  near  -6  MHz  and  indicate  a 
radial  wind  velocity  of  +6  m/sec  (positive  velocity  away  from  the  lidar).  Note  the 
interesting  velocity  structure  within  the  first  2.5  km,  especially  the  5  MHz  (5  m/sec) 
shear  from  1 .4  km  to  1 .9  km. 

In  Figure  1.7b,  the  first  Doppler  spectral  moment  (velocity)  is  plotted  versus  range  in  a 
rectangular,  side-looking  PPI  display.  For  this  plot,  the  lidar  scanned  through  90°  in 
azimuth  with  the  elevation  angle  held  constant  at  +3°  (this  data  was  collected  on  a 
different  day  from  that  in  Figure  1.5).  The  color  lookup  table  at  the  bottom  of  the 
screen  extends  from  -8  m/sec  (blue)  to  +8  m/sec  (red).  (Positive  velocities  are  defined 
as  those  away  from  the  lidar).  Over  the  90°  azimuth  scan,  the  radial  wind  velocity  shifts 
from  values  near  +3  m/sec  to  values  near  -7  m/sec.  From  an  azimuth  angle  of  90°  to 
an  azimuth  angle  of  135°,  interesting  small-scale  (1  to  2  m/sec  over  ranges  of  200  to 
400  m)  wind  eddies  are  clearly  present. 

Figure  1.9  shows  the  second  set  of  data  which  was  taken  with  the  one  micron  system 
set  to  scan  a  sector  in  a  vertical  plane.  This  system  operates  at  a  higher  prf  (10  Hz) 
with  80  millijoule  pulse  energy  and  used  a  0.2  microsecond  pulse.  This  system  is 
beginning  to  approach  the  36  Hz  system  requirement  quoted  above.  In  Figure  1.9  the 
data  have  been  processed  to  show  both  Doppler  mean  velocity  and  backscatter 
intensity  vs  range  for  successive  pulses  as  the  scan  executed  a  nodding  pattern  in  the 
vertical  direction  The  mean  velocity  was  calculated  using  a  pulse  pair  algorithm  and 
each  point  represents  the  average  of  about  10  pulses.  The  scan  cycle  is  clearly 
evident  in  both  the  intensity  and  the  velocity  data.  (The  strong  stationary  feature 
located  at  about  1  km  range  is  an  artifact  echo  which  resulted  from  a  system  fault  at  the 
time  of  the  data  collection). 

These  data  were  taken  during  an  experiment  to  measure  aircraft  wake  velocities. 
Figure  1.10  shows  a  contour  plot  of  the  line  of  sight  velocity  seen  in  the  vertical  plane 
with  this  system  (a  different  data  set  than  that  in  Figure  1.9).  The  quantity  shown  is  the 


14 


16 


;310  AZIM  14-4.057  ELEV  003. 0C9 


velocity  -  PP 


Intensity 


Figure  1.9  Pulse  -  pair  estimates  of  velocity  and  backscattered  intensity  vs 

range  at  1  micron.  Every  tenth  pulse  is  shown  as  the  elevation  angle  of 
the  look  direction  is  scanned  through  a  9  degree  sector  between  18  and 

27  degrees.  The  pulse  covariance  is  averaged  over  10  successive 

pulses  (1  second  at  the  10  Hz  prf)  before  application  of  the  pulse  pair 
alogorithm.  The  range  extent  shown  is  approximately  1.5  km.  The  zigzag 
pattern  is  due  to  the  oscillating  scan  (20  second  period-linear  profile)  as 
the  line  of  sight  sweeps  through  the  vertically  sheared  atmosphere  The 
strong  feature  near  1  km  is  an  artifact. 
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first  moment  of  the  Doppler  spectrum  in  successive  30  meter  range  gates  and  is  plotted 
as  a  function  of  range  and  crossrange  in  the  measurement  plane.  The  strong  feature 
evident  in  the  lower  right  of  the  plot  is  the  vortex  wake  of  a  Boeing  707  aircraft  which 
passed  through  the  scan  plane  a  few  seconds  earlier.  It  is  evident  from  this  figure  that 
velocity  structures  having  scales  as  small  as  25  meters  can  be  resolved  with  this  one 
micron  system. 

The  preliminary  measurements  described  above  demonstrate  that  both  the  two  micron 
and  the  one  micron  system  have  the  basic  capability  to  acquire  velocity  data  with  the 
characteristics  needed  for  the  boundary  layer  analysis.  The  principal  additional  need 
for  a  full  PBL  analysis  is  to  achieve  pulse  repetition  rates  up  to  several  hundred  Hz. 
Diode  pumped  laser  systems  are  currently  being  developed  at  CTI  that  have  this  latter 
capability. 


System  Performance  and  SNR  Model  Validation 

The  accuracy  with  which  the  wind  velocities  can  be  inferred  from  lidar  measurements 
depends  primarily  on  two  quantities:  the  system  SNR  and  the  Doppler  spectrum  width. 
Both  of  these  quantities  depend  on  the  atmospheric  scattering  and  propagation 
properties  as  well  as  on  the  lidar  system  parameters.  CTI  maintains  an  extensive  set  of 
simulation  and  prediction  codes  for  predicting  the  performance  of  coherent  lidar 
systems  operating  against  atmospheric  targets.  The  utility  of  these  models  relies 
heavily  on  how  well  atmospheric  parameters  are  known  both  horizontally  and  vertically. 
The  atmospheric  parameters  of  interest  are  the  volume  backscatter  coefficient,  the 
atmospheric  extinction  coefficient,  and  the  refractive  turbulence  structure  coefficient. 

Table  1  shows  an  example  of  a  comparison  of  measured  values  of  SNR  with  those 
predicted  by  the  model.  The  2.09  mm,  21  mj,  5  Hz,  10  cm  aperture  (8  cm  e*2  intensity 
diameter)  system  was  employed  for  these  measurements  during  a  recent  Wright 
Laboratory-sponsored  program.  The  skies  were  overcast  during  the  measurements,  so 
that  refractive  turbulence  effects  should  be  small  (although  no  direct,  quantitative 
measurement  of  refractive  turbulence  was  made).  For  comparison  purposes,  model 
predictions  for  the  standard  refractive  turbulence  structure  coefficient  model  (Hufnagel) 
are  given  in  addition  to  this  same  model  reduced  by  a  factor  of  0.01.  As  expected,  the 
latter,  'weak'  turbulence  model  predictions  better  match  the  measured  values.  The 
agreement  is  quite  good  for  all  elevation  angles  and  all  ranges. 


1.6  DEVELOPMENT  PLANS 

The  demonstration  measurements  acquired  during  this  study  together  with  the 
predictions  of  performance  and  the  development  of  the  vector  wind  and  water  vapor 
correlation  processing  indicate  that  three  dimensional  profiling  of  the  vector  wind  and 
water  vapor  concentration  from  a  single  lidar  station  of  the  time  dependent  planetary 
boundary  layer  at  resolutions  of  25  to  50  meters  should  be  possible.  Two  new 
processing  concepts  have  been  proposed  for  this  analysis  and  these  will  need  to  be 
evaluated.  A  program  to  develop  such  a  boundary  layer  probe  is  outlined  in  the  last 
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Table  1 

SNR  MODEL  VALIDATION  —  PRELIMINARY 
OVERCAST  SKIES 


ELEVATION 

ANGLE 

RANGE  (km) 

PRED.  SNRn  (dB) 
Standard  Cn2 

PRED.  SNRn  (dB) 
0.01 ‘(Stand.  Cn2) 

MEAS.  SNRn  (dB) 

0° 

1 

21.0 

27.5 

25.1 

2 

12.7 

25.2 

21.5 

4 

3.6 

20.2 

17.1 

6 

-2.0 

16.0 

14.7 

8 

-6.2 

12.6 

12.2 

2° 

1 

25.1 

27.4 

28.1 

2 

19.6 

25.0 

23.2 

4 

11.9 

20.1 

17.5 

6 

6.4 

15.9 

15.7 

8 

1.9 

12.3 

13.8 

5° 

1 

25.8 

27.0 

27.2 

2 

21.0 

24.4 

23.6 

4 

12.7 

18.8 

19.3 

6 

6.3 

14.0 

14.3 

8 

1.1 

9.9 

10.5 

10° 

1 

25.7 

26.4 

27.9 

2 

20.8 

23.2 

25.4 

4 

11.5 

16.8 

11.9 

6 

4.1 

11.2 

11.0 

8 

-1.8 

6.5 

10.7 

20° 

1 

24.8 

25.3 

28.7 

2 

19.2 

21.2 

13.9 

4 

8.3 

13.2 

11.4 

90° 

2 

12.7 

14.4 

14.7 

\  =  2.09  urn;  D  =  8  cm  (e-2  intensity  diameter);  Ej  =  21  mJ;  x  =  220  ns; 

T|  =  0.1 
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llijoule,  0.3  microsecond 
lere,  light  turbulence,  0.1 


section  of  the  report.  The  program  requires  the  construction  or  acquisition  of  a  1  or  2 
watt  mean  power,  pulsed  diode  pumped  laser  that  can  operate  at  several  hundred  Hz. 
The  expected  SNR  for  a  candidate  boundary  layer  measuring  system  (approximately  7 
millijoules  at  2.1  microns  with  a  0.3  microsecond  pulse  and  a  0.1  meter  aperture)  is 
shown  in  Figure  1.8  as  a  function  of  range  and  altitude(midlatitude  winter  clear 
atmosphere,  light  turbulence).  SNR  values  exceed  10  db  to  ranges  of  4  km  and 
altitudes  exceeding  1  km  and  provide  an  accurate  velocity  measurement  capability  in 
this  region. 


2.0  MEASUREMENT  CONCEPTS 
2.1  OBJECTIVES 

The  objective  of  the  study  has  been  to  evaluate  the  capability  for  a  scanning  coherent 
pulsed  laser  to  acquire  high  resolution  mapping  of  the  three  dimensional  structure  of 
wind  velocity  and  water  vapor  in  the  turbulent  planetary  boundary  layer.  This 
information  is  needed  to  obtain  a  better  understanding  and  modeling  of  the  turbulent 
transport  mechanisms  within  the  boundary  layer  to  support  the  development  and 
improvement  of  subgrid  models  for  the  transport  that  can  be  used  in  large  scale  eddy 
simulations  of  the  turbulent  motion. 


2.1.1  MEASUREMENT  CONCEPTS 

A  coherent  lidar  senses  the  component  of  the  motion  of  the  atmosphere  parallel  to  its 
instantaneous  look  direction.  To  infer  vector  winds  it  is  necessary  either  to  have 
multiple  looks  at  each  point  in  space  from  different  directions,  or  to  invoke  prior 
knowledge  that  the  atmospheric  motion  is  sufficiently  constrained  b>  the  requirements 
of  mass,  momentum  and  energy  conservation  that  the  full  vector  field  can  be 
determined  from  measurements  of  a  single  component.  We  assume  that  a  3D  volume 
of  atmosphere  is  scanned  from  a  ground  station  or  an  aircraft.  To  determine  the  vector 
wind  at  each  point  in  space  without  making  any  assumptions  about  the  atmosphere, 
three  separated  lidars  would  be  needed  to  deduce  the  3D  wind  vector  at  each 
observation  point.  It  is  the  conclusion  of  the  present  study  that  winds  in  the  actual 
atmosphere  are  sufficiently  constrained  by  the  rules  of  fluid  mechanics  that  the  same 
capability  can  be  achieved  under  most  circumstances  with  a  single  lidar  which  senses 
only  the  instantaneous  component  of  velocity  parallel  to  the  line  of  sight.  By  making 
such  measurements  quasi  continuously  in  time  over  a  3D  volume  of  space  and  by 
making  use  of  known  or  assumed  constraints  on  the  motion,  the  vector  information  can 
be  retrieved. 

This  basic  approach  has  been  well  developed  in  radar  meteorology  where  data 
acquired  with  two  separated  radars,  viewing  a  common  storm  volume  is  combined  with 
the  constraint  of  fluid  continuity  to  yield  the  3D  distribution  of  wind  vectors  throughout 
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the  storm  environment  ^  Similar  applications  of  momentum  constraints  have  been 
applied  to  radar  interrogation  of  storm  clouds  from  fast  moving  aircraft(3,4).  Based  on 
the  results  of  the  present  study,  we  feel  that  extensions  of  these  techniques  are 
appropriate  for  lidar  observations  of  the  planetary  boundary  layer.  In  particular  we  feel 
that  3  dimensional  scanning  measurements  made  from  a  single  platfrom,  either  ground 
or  airborne,  can  be  processed  to  yield  a  3  dimensional  description  of  the  wind  vector 
field  in  the  turbulent  boundary  layer.  The  ability  to  do  this  relies  on  having  three 
elements  in  the  processing:  1)  an  accurate  model  of  the  measurement  process  (i.e., 
how  the  speckled  Doppler  data  is  related  to  the  scattering  and  motion  of  the  aerosol 
scatterers  in  the  focal  volume),  2)  a  realistic  mathematical  description  of  the  three 
dimensional  fluid  dynamics  of  the  hydrodynamic  constraint,  and  3)  a  well  founded 
algorithm  for  merging  these  models  with  the  data. 

In  this  report  we  will  outline  the  proposed  measurement  and  processing  concepts  and 
provide  estimates  and  simulations  of  expected  performance  for  candidate  measuring 
geometries. 

The  correlation  of  the  vector  wind  field  with  the  distribution  of  passive  constituents, 
particularly  water  vapor,  is  key  to  the  understanding  of  energy  transport  in  the 
convectively  turbulent  boundary  layer.  DIAL  (Differential  Absorption  Lidar)  is  a 
technique  which  senses  spectral  differences  due  to  water  vapor  absorption  in  displays 
of  backscattered  lidar  energy.  In  principle  data  from  a  scanning  DIAL  system  can  be 
processed  to  yield  a  description  of  the  3D  spatial  structure  of  the  water  vapor  mixing 
ratio  at  the  same  scales  as  the  wind  field  measurements.  However,  DIAL  as  it  is 
normally  implemented,  is  very  insensitive  to  structure  scales  that  are  much  less  the 
overall  scale  of  the  absorbing  path.  Thus  it  is  most  useful  for  describing  the  large  scale 
distributions  of  the  absorbing  species.  In  this  study  we  have  developed  a  modified 
DIAL  concept  which  relies  on  correlations  of  the  water  vapor  fluctuations  with 
fluctuations  in  the  density  of  natural  aerosols  over  a  range  of  spatial  scales.  By  using 
DIAL  to  sense  the  correlation  between  aerosol  and  water  vapor  fluctuations  at  large 
scales  plus  direct  measurements  of  aerosol  fluctuations  at  the  smaller  turbulent  scales, 
the  water  vapor  fluctuations  at  both  large  and  small  scales  can  be  estimated.  When 
combined  with  simultaneous,  colocated  velocity  measurements,  the  turbulent  fluxes  of 
moisture  and  latent  heat  can  be  evaluated. 

In  this  report  we  will  identify  a  suitable  measurement  and  processing  system  for 
measuring  the  combined  water  vapor  -  velocity  field  and  outline  a  candidate 
measurement  program  to  confirm  the  concepts. 


2.1.2  MEASUREMENT  OF  VECTOR  WINDS 

We  consider  two  basically  different  measurement  concepts:  one  from  an  aircraft  that 
flies  over  or  under  the  region  to  be  surveyed,  and  the  other  a  ground  installation  that 
views  the  atmosphere  above  and  around  the  lidar  installation.  To  deduce  a  three 
component  wind  vector  at  each  point  in  space  from  measurements  of  one  component, 
at  least  two  additional  conditions  must  be  imposed.  In  the  present  analysis  these  will 
be  the  requirements  of  continuity  and  of  momentum  conservation. 
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Mass  conservation  requires  that  the  deduced  velocity  vector  field  be  divergence  free, 
i.e.,  that  the  equation  of  continuity  is  satisfied: 


dii  +  dv  dw 
dx  dy  dz 


(1) 


where  u,v,w  are  the  three  velocity  components  corresponding  to  the  directions  x,y,z  (we 
use  the  convention  that  z  and  w  are  in  the  vertical  direction).  Measurements  of  two 
components  of  velocity  over  the  entire  3D  space  plus  this  equation  allow  the  full 
reconstruction  of  the  vector  field. 


This  method  has  been  used  with  great  success  in  radar  measurements  of  storm 
structure  (Dual  Doppler  method).  Here  measurements  of  the  line  of  sight  velocity  of  a 
storm  region  are  made  simultaneously  from  two  different  ground  stations.  When 
combined  with  Eq.  (1)  together  with  the  boundary  condition  that  the  vertical  velocity 
vanish  at  the  ground  surface,  the  3D  distribution  of  vector  winds  can  be  reconstructed. 

An  example  of  such  a  3D  vector  wind  retrieval  is  shown  in  Figure  2.1.  Since  the  radars 
view  the  atmosphere  continuously,  a  3D  motion  picture  of  the  vector  wind  field  can  be 
constructed  by  this  method.  This  technique  has  been  very  useful  in  eliciting 
understanding  of  the  detailed  dynamics  of  severe  storms  and  their  generation  of  low 
altitude  wind  shear.  In  principle  the  same  procedure  can  be  applied  to  lidar  using  two 
spatially  separate  systems  that  conduct  a  coordinated  scan  of  a  common  atmospheric 
volume. 


The  Dual  doppler  method  requires  two  separate  ground  stations.  An  approximation  to 
such  a  multi  view  system  with  a  single  radar  or  lidar  can  be  accomplished  on  board  a 
fast  moving  aircraft.  This  technique  has  been  proposed  and  implemented  by  a  number 
of  experimenters  for  radar  observation  of  local  storms  from  aircraft. 

In  this  flyby  mode  the  aircraft  scans  the  3D  field  above  or  below  (or  around)  the  aircraft 
as  it  flies  past  the  atmospheric  volume  to  be  scanned.  As  long  as  the  aircraft  speed  is 
much  greater  than  the  atmospheric  motions  being  observed,  the  atmosphere  can  be 
assumed  frozen  in  its  motion  for  the  period  required  for  the  aircraft  viewing  aspect  to 
change  substantially.  Roughly  speaking,  the  aircraft  speed  should  be  many  times 
larger  than  the  product  of  the  rms  atmospheric  velocity  fluctuations  and  the  ratio  of  the 
overall  dimension  of  the  absorbing  region  to  the  spatial  scale  of  the  fluctuations  of 
interest: 


V»(R/l)v 


(2) 
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JO  M  JO  -  <acl  tel. 

Vortlcol  ond  Horlioncol  Wli*4»,  Ull  CJT  OoMociiviiy.  HII  CST 

The  latest  instrument  in  severe  thunderstorm  research  ts  dual  Doppler  radar.  Data  from  two  Doppler  radar 
sets  aimed  at  a  single  thunderstorm  can  be  used  to  reconstruct  the  three-dimensional  (low  patterns  as  shown  (or  a 
storm  that  occurred  in  Oklahoma  on  June  8.  1974  These  data  revealed  the  development  o!  a  double  vortex  inside 
a  supercell  thunderstorm  as  a  tornado  extended  lo  the  ground  The  radar  reflectivity  patterns  at  2  km  show  the  hook 
echo  m  relationship  to  the  rotation  of  the  mesocyclomc  vortex  that  contributes  to  its  shape  (from  J  R  Eagieman  and 
W.  C.  bn.  "Severe  Thunderstorm  Internal  Structure  from  Dual  Doppler  Radar. '  J  Appl.  Meteorology.  October  1977) 


Figure  2.1  Dual  Doppler  Radar  Reconstruction 
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Here  we  take  the  range  R  as  a  measure  of  the  outer  scale  of  the  viewed  region  and  / 
and  v  to  be  the  fluctuation  scale  and  velocity  respectively.  This  condition  states  that,  in 
order  to  freeze  the  motion,  the  rotation  rate  of  the  viewing  line  of  sight  (V/R)  must 
greatly  exceed  the  physical  rotation  rate  of  typical  fluid  elements  at  the  measurement 
scale  (v/l).  For  V=150  m/s  and  v=5  m/sec,  this  requirement  will  not  be  met  for 
fluctuation  scales  much  smaller  than  a  third  of  the  range.  Thus  application  of  this 
method  must  be  viewed  with  caution  for  following  intensely  turbulent  features  and  may 
in  fact  not  be  very  reliable  for  most  turbulent  structures.  It  is  most  appropriate  for 
diagnosing  very  weak  turbulence  structures  (<1  m/s  fluctuations)  or  for  eliciting  the 
large  scale  motions  of  storm  regions. 

From  a  single  ground  station  the  same  'frozen'  structure  approach  can  be  applied  if  it  is 
assumed  that  the  atmosphere  drifts  past  the  observing  station  without  sensible  change 
in  its  small  scale  structure.  Here  the  mean  drift  velocity  replaces  the  aircraft  velocity  in 
Eq.  (1).  The  region  of  applicabilty  of  this  method  is  very  limited  for  a  ground  station 
unless  it  is  known  a  priori  that  small  scale,  long-lived  advecting  patterns  are  present. 

A  better  approach  to  the  interpretation  of  data  from  both  a  single  ground  station  or  from 
an  aircraft  platform  can  be  achieved  if  a  more  physical  account  of  the  advection 
process  is  used.  Here,  instead  of  invoking  a  'frozen'  turbulence  assumption,  a 
requirement  that  the  observations  be  consistent  with  conservation  of  momentum  is  to 
be  imposed.  This  constraint  can  be  imposed  by  demanding  that,  in  addition  to 
satisfying  the  continuity  equation,  the  inferred  velocity  field  also  satisfy  momentum 
conservation  equations.  In  the  limit  of  rapid  but  steady  advection  of  a  slowly  changing 
atmosphere  this  constraint  reduces  to  the  frozen  advection  assumption.  However, 
since  it  is  applicable  to  the  atmosphere  as  observed  from  platforms  moving  at  any 
speed,  including  zero,  it  has  a  much  greater  range  of  applicability.  This  approach  has 
been  used  in  this  study.  The  mathematics  of  its  invocation  are  discussed  in  detail  in 
Section  3.1. 


2.1.3  MEASUREMENT  OF  WATER  VAPOR  USING  DIAL 

Transport  of  water  vapor  in  the  atmosphere  is  a  vital  component  of  the  atmospheric 
dynamics.  Latent  heat  release  or  extraction  as  water  condenses  or  evaporates  is  an 
important  and  often  dominating  source  of  temperature  change  and  the  buoyant  forces 
that  drive  the  convective  motions  in  unstable  boundary  layers.  The  upward  transport  of 
energy  is  determined  by  correlations  between  water  vapor  fluctuations  and  fluctuations 
in  the  vertical  velocity.  The  goal  of  the  present  study  is  to  evaluate  the  potential  of 
using  lidar  to  monitor  simultaneously  both  of  these  quantities  at  scales  as  small  as  20 
meters.2  Differential  absorption  of  the  lidar  energy  scattered  from  natural  aerosols 
(DIAL)  is  a  technique  to  deduce  the  presence  and  spatial  distribution  of  an  absorbing 
species  by  comparing  the  amplitude  of  the  backscattered  light  at  two  different,  but 
closely  spaced,  wavelengths  where  the  extinction  coefficients  differ  substantially  but  for 
which  the  scattering  coefficients  are  virtually  identical. 


2T1ic  estimates  presented  later  in  this  report  indicate  that  estimation  of  the  velocity  with  submeter  per 
second  precision  may  limit  the  spatial  resolution  to  somewhat  larger  values  (25  to  50  meters). 


26 


Figure  1.4a  shows  the  absorption  coefficient  of  a  midlatitude  summer  atmosphere  at 
ground  level  in  the  neighborhood  of  2.01  microns;  Figure  1.4b  shows  the  absorption 
coefficient  due  just  to  water  vapor  alone.  With  a  Tm:YAG  laser  any  wavelength  in  the 
spectral  region  indicated  can  be  selected  by  electronic  tuning.  In  the  DIAL  technique 
using  a  pulsed  laser,  two  wavelengths  are  alternately  selected  at  successive  pulses  or 
intervals:  one  in  the  vicinity  of  the  center  of  an  absorbing  line  and  one  immediately 
adjacent  where  the  absorption  coefficient  is  substantially  smaller.  To  achieve  maximum 
sensitivity,  the  absorption  coefficient  should  be  chosen  so  that  the  extinction  in  the 
absorbing  channel  is  about  a  factor  2  or  3  at  the  center  of  the  region  to  be  diagnosed; 
i.e.,  so  that  the  absorption  coefficient  is  comparable  to  the  inverse  of  the  mean  range. 

As  an  example,  one  such  line  for  the  Tm.YAG  laser  is  marked  in  Figure  1.4.  The 
center  of  the  line  is  located  at  2.016974  microns.  Since  the  line  half  width  at  half 
maximum  is  2.8  GHz/atm,  a  separation  greater  than  3  GHz  is  suitable  for  the  non¬ 
absorbing  wavelength.  At  2.017713  microns  (5.8  GHz  separation)  there  is  a  minimum 
in  the  total  absorption  equal  to  about  0.1  per  km.  At  one  atmosphere  pressure,  the 
absorption  coefficent  at  the  line  center  is  about  3.4  per  kilometer  (midlatitude  summer 
atmosphere).  Thus  two  wavelengths  can  always  be  selected  in  this  region  that  will  be 
optimal  (  i.e.,  for  which  the  absorption  coefficient  of  the  stronger  line  is  of  order 
1/range)  for  a  DIAL  measurement  of  water  vapor  for  any  nominal  range  between  a  few 
hundred  meters  to  more  than  5  km. 

Two  pulses,  closely  spaced  in  time,  are  to  be  transmitted,  one  at  each  wavelength.  In 
expectation  the  return  signal  intensities  at  the  two  wavelengths  will  have  the  form 


A00  =  A(*)0(*)exp(-2j*,(x>k) 

0 

(3) 

X 

AO)  =  A  0)00)  exp  (-2  J  k2(x)dx) 

o 


Here  700)  contains  the  geometry  factors  associated  with  the  illumination  and  reception 
and  P(x)  is  the  local  backscatter  coefficient  (both  assumed  identical  at  the  two 
wavelengths).  The  local  absorption  coefficient  difference  (and  therefore  the  local  water 
vapor  density  or  mixing  ratio)  is  evaluated  as  a  function  of  range  (in  the  absence  of 
noise)  by  differentiating  the  logarithm  of  the  ratio  of  these  two  signals: 


Ak(x)  =  ki  (x)  -  k2  (x)  =  ^  i w(/2  (x)  /  /,  (x)) 

2  ax 


(4) 
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With  finite  noise  a  least  squares  formula  should  be  used: 


A k(x)=)-^-ln(<  I2(x)Ii(x)>  /  </,(x)3  >) 
2  ax 


(5) 


In  Eq.(5)  the  brackets  (<>)  represent  an  average  over  successive  samples  or 
neighboring  range  gates. 

The  DIAL  measurement  is  limited  by  the  noise  in  the  intensity  channel,  especially  when 
high  spatial  resolution  is  desired.  This  noise  is  due  to  speckle  at  high  SNR  and  is 
photon  shot  noise  at  low  SNR.  The  expected  fractional  error  (standard  deviation)  of  the 
measured  absorption  coefficient  is  given  approximately  by  (5) 


5k  _  1  Jl  +  SNR2 

k  ~  kAyfN  SNR 


Here  A  is  the  resolution  width  and  N  the  number  of  independent  samples  contributing  to 
the  estimate.  At  large  SNR  and  for  an  optimum  level  of  absorption(k~1/L),  the  rms 
precision  (  dk )  of  the  inferred  absorption  coefficient  is  given  by 


(L/A) 
(8k  Ik)2 


For  meteorological  estimates  a  precision  of  the  estimate  of  the  local  water  vapor 
content  of  the  order  of  10%  (standard  deviation)  or  better  is  probably  necessary.  Using 
DIAL,  even  with  optimum  choice  of  absorption  coefficient,  a  very  large  number  (104  to 
10^)  of  independent  samples  per  resolution  element  will  be  needed  to  achieve  a  useful 
measurement  of  water  vapor  fluctuations  at  scales  less  than  a  tenth  of  the  range.  Thus 
the  DIAL  technique  is  best  applied  to  estimating  relatively  coarse  grained  averages. 

It  is  very  difficult  to  achieve  such  high  levels  of  averaging  while  maintaining  high 
detection  efficiency  with  a  heterodyne  detection  system.  Without  defocusing  each 
pulse  is  limited  to  providing  one  or  at  most  a  few  independent  samples  (depending  to 
the  degree  velocity  spread  within  the  pulse  length).  In  order  to  achieve  such  high  levels 
of  averaging  it  is  necessary  to  detect  a  very  large  number  of  incoherent  scatterers 
within  each  resolution  field  to  reduce  the  speckle  noise.  Illumination  with  spectrally  and 
angularly  defocused  or  broad  laser  pulses  and  using  direct  (non-coherent)  detection  of 
the  total  backscattered  energy  is  required.  Even  so  the  averaging  requirements  can  be 
severe  when  high  spatial  resolution  is  desired  (see  Section  2.1.4). 

In  view  of  the  above  considerations  we  have  concluded  that  the  DIAL  technique  as 
conventionally  implemented  does  not  provide  a  useful  method  with  a  coherent  lidar 
using  heterodyne  detection  for  measuring  spatial  distributions  of  water  vapor 
fluctuations  at  small  turbulence  scales.  DIAL  is  most  appropriate  to  the  sensing  of 
moderate  to  large  scale  features  (>10%  of  the  overall  dimension). 
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Alternative  strategies  include  direct  detection  of  the  backscattered  energy  from  non 
coherently  illuminated  scatterers  in  each  resolution  volume  or  correlation  of  the  water 
vapor  fluctuations  with  the  more  readily  detected  aerosol  fluctuations.  These  will  be 
discussed  next. 


2.1.4  MEASUREMENT  PRECISIONS  :  Aerosol  scattering  intensity  and  DIAL 

It  is  informative  to  examine  how  the  energy  and  prf  requirements  for  a  DIAL 
measurement  of  the  density  distribution  of  an  absorbing  species  and  direct 
measurements  of  the  local  backscattering  coefficient  depend  on  the  spatial  resolution 
of  the  measurement. 

We  assume  that  a  3D  volume  of  dimension  L  is  to  be  scanned  in  a  specified  time  T  with 
a  3D  resolution  A.  Two  sources  of  noise  contribute:  speckle  and  photon  shot  noise. 

Photon  Shot  Noise  Only 

The  total  number  of  photons  that  must  be  incident  on  the  interrogated  volume  during 
the  dwell  time  to  achieve  the  desired  rms  detection  precision  e  of  the  scattering 

coefficient  is  proportional  to 


N 


scatt 


(8) 


where  G  is  a  factor  determined  by  the  collection  geometry. 

In  a  ideal  DIAL  measurment  we  would  choose  the  absorption  coefficient  to  be 
approximately  equal  to  1/L.  For  this  choice  the  number  of  photons  that  need  to  be 
incident  to  achieve  the  same  precision  in  the  measurement  of  absorption  coefficient  is 
equal  to 


1  1 


e2  M 


(9) 


i.e.  a  factor  of  (L/A)2  larger.  This  additional  factor  in  DIAL  results  from  the  fact  that  the 
DIAL  measurement  must  difference  adjacent  range  cells  to  get  the  absorption  signature 
whereas  the  backscatter  is  a  direct  measurement.  Thus  when  photon  noise  dominates 
it  is  much  easiest  to  measure  the  mean  scattering  coefficient  than  it  is  to  measure  the 
absorption  coefficient.  When  speckle  dominates,  this  proportionality  and  these 
formulae  still  hold  as  long  as  N  is  interpreted  as  being  the  number  of  independent 
samples  acquired,  not  the  number  of  photons  incident. 


Because  of  the  very  strong  dependence  of  required  power  and  number  of  independent 
samples  on  spatial  resolution  for  the  DIAL  measurement  (iV*l/A5)  we  can  expect 
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that  DIAL  measurements  will  be  restricted  to  relatively  coarsed  grained  averages  in 
space. 


2.1.5  DIRECT  DETECTION  DIAL 

Direct  (incoherent)  detection  lidar  systems  detect  the  backscattered  electric  field 
intensity  and  do  not  retain  phase  information.  Because  they  are  intensity-only 
detection  systems,  wavefront  distortions  due  to  atmospheric  refractive  turbulence  only 
negligibly  reduce  system  SNR.  Consequently,  these  systems  are  able  to  employ 
significantly  larger  apertures  (50  cm  and  more)  without  turbulence-induced  coherence 
loss  effects.  Furthermore,  by  using  a  transmit  aperture  that  is  much  smaller  than  the 
receiver  aperture,  direct  detection  systems  are  able  to  reduce  and  effectively  eliminate 
laser  speckle  effects.  The  transverse  size  of  laser  speckle  lobes  at  the  lidar  is 
approximately  equal  to  the  transmit  aperture  size.  The  larger  receiver  aperture  collects 
backscattered  radiation  from  several  speckle  lobes  simultaneously  and  thereby 
'aperture  averages'  out  the  speckle-induced  intensity  fluctuations.  The  number  of 
speckle  lobes  that  can  be  averaged  with  direct  detection  lidar  systems,  or  the  degree  M 
of  aperture  averaging  available,  is  proportional  to  the  ratio  of  the  receiver  aperture  area 
to  the  transmit  aperture  area.  For  example,  a  3  cm  transmit  aperture  with  a  48  cm 
receiver  aperture  achieves  M  ~  256.  To  achieve  this  same  sort  of  speckle  reduction,  a 
heterodyne  detection  system  would  need  to  average  256  pulses  or  equivalently, 
employ  a  256  element  detector  array. 

Water  vapor  measurements  with  direct  detection  DIAL  has  been  performed  most 
notably  at  ruby  (694  nm),  CO2  (10  pm),  and  Alexandrite  (725-730  nm)  laser 
wavelengths.  One  of  the  primary  difficulties  with  these  all  of  these  systems  is  that  they 
cannot  be  made  sufficiently  compact  while  at  the  same  time  providing  accurate,  high 
spatial  and  range  resolution,  water  vapor  measurements  over  sufficiently  large  areas. 


2.1.6  CORRELATION  DIAL 

An  alternate  strategy  for  estimating  the  water  vapor  fluctuations  utilizes  the  expected 
space  time  correlation  that  should  exist  in  a  turbulent  flow  between  ail  passively 
advected  constituents,  in  particular  water  vapor  and  ambient  aerosols.  In  this  technique 
DIAL  measurements  at  large  scales  are  used  to  interpret  aerosol  density 
measurements  at  small  scales  in  terms  of  estimates  of  water  vapor  fluctuations. 

The  basic  assumption  is  that  all  passive  constituents  that  are  advected  in  the  turbulent 
atmosphere  by  the  same  flow  mechanics  will  have  a  similar  and  coherent  distribution  of 
fluctuations  over  a  range  of  scales.  This  requirement  probably  limits  the  applicability  of 
the  method  to  clear  or  hazy  atmospheres  that  are  free  from  clouds  or  fog.  Evaporation 
and  condensation  of  water  droplets  near  or  in  clouds  may  degrade  the  correlation 
between  the  net  backscattering  coefficient  and  the  water  vapor  absorptivity. 
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Given  that  the  uniform  passive  advection  assumption  is  valid,  a  transfer  function 
between  water  vapor  fluctuations  and  aerosol  density  fluctuations  can  be  established  at 
large  scales  where  fluctuations  of  both  constituents  can  be  measured  with  reasonable 
accuracy.  As  was  shown  in  Section  2.1.4  the  fluctuations  of  the  scattering  intensity  of 
the  aerosols  are  much  more  readily  measured  at  small  scales  with  lidar  than  are  the 
absorption  fluctuations.  By  assuming  that  the  transfer  function  between  these  two 
quantities  is  independent  of  scale,  the  small  scale  water  vapor  structure  can  be 
deduced  from  measurements  of  the  aerosol  structure  by  multiplying  by  the  transfer 
function  derived  from  measurements  at  large  scale. 

For  a  coherent  lidar,  the  measurements  of  fluctuations  of  aerosol  density  at  high  SNR 
obey  Rayleigh  statistics  due  to  speckle.  For  a  quiescent  atmosphere,  the  correlation 
length(&)  is  the  pulse  or  resolution  length(A).  For  a  turbulent  atmosphere  the 
correlation  length  may  be  smaller  and  is  determined  by  the  velocity  spread  within  the 
pulse  length. 

The  expected  fractional  standard  deviation  of  the  local  scattering  coefficient  varies  as 


5p_  [&  1  Jl  +  SNR2 
p  ~  V  A  yfN  SNR 


The  number  cf  samples  required  to  achieve  a  given  precision  in  the  estimate  of  the 
scattering  coefficient  (and  therefore  of  the  water  vapor  density)  is  given  by 


WP)2 


(11) 


Since  &  <  A  in  general,  the  number  of  samples  required  by  this  method  will  be  much 

less  than  that  for  direct  DIAL  (by  the  factor  (LI 5s)2).  This  should  allow  estimation  of 
water  vapor  fluctuations  down  to  scales  limited  only  by  a  requirement  to  achieve  an 
adequate  single  range  gate  SNR. 


2.2  APPROACH 

The  proposed  approach  for  measuring  water  vapor  fluctuations  simultaneously  with 
wind  vector  fluctuations  may  be  summarized  as  follows. 

-  Radial  winds  and  backscattering  coefficients  are  measured  over  a 
3D  volume  as  a  function  of  time  using  a  single  scanning  coherent 
lidar. 

-  Constraints  of  continuity  and  momentum  conservation  are  used  to 
infer  the  distribution  of  the  vector  wind  over  this  volume. 


31 


-  Using  one  of  several  spatial  horizontally  averaging  schemes,  mean 
water  vapor  and  backscattering  coefficients  are  measured  and 
correlated  to  yield  three  vertical  profiles: a  height  dependent  transfer 
function,  the  horizontally  averaged  mean  water  vapor  mixing  ratio, 
and  the  horizontally  averaged  mean  backscattering  coefficient.  These 
measurements  are  to  be  carried  out  before,  during,  and/or  after  the 
wind  vector  measurements. 

-  The  water  vapor  -  aerosol  transfer  function  is  used  to  convert  the 
scattering  coefficient  measurements  at  small  spatial  scales  to  a  3D 
distribution  of  water  vapor  concentration. 


3.0  PERFORMANCE  MODELING 

3.1  SIGNAL  PROCESSING  AND  VECTOR  WIND  RETRIEVAL 

3.1.1  KALMAN  FILTER  RETRIEVAL  PROCESSING  THEORY 

A  full  3D  approximation  of  the  Kalman  filter  wind  field  retrieval  analysis  has  been 
formulated  analytically(  see  Appendix).  In  the  following  section  3.1.2  we  summarize  the 
mathematical  formalism  for  retrieval  of  vector  winds  for  single  station  Doppler 
measurements  of  the  line  of  sight  velocity. 

3.1.2  VECTOR  WINDS  FROM  DOPPLER  DATA 

Two  formulations  have  been  analysed:  an  Eulerian  form  in  which  velocity  and  vorticity 
are  evaluated  and  tracked  on  a  fixed,  prespecified  grid,  and  a  Lagrangian  form  in 
which  a  cloud  of  discrete  vortex  elements  or  threads  are  followed.  Flow  charts  of  these 
procedures  are  shown  in  Figures  3.1  and  3.2.  The  specific  analysis  for  incorporating 
the  lidar  observations  into  the  fluid  dynamic  equations  describing  the  momentum  and 
mass  conservation  equations  of  the  atmospheric  boundary  layer  follows  previous 
concepts  but  is  fundamentally  a  new  formulation. 

The  formulation  for  the  Lagrangian  mode,  whose  basic  equations  are  described  in  this 
section,  has  a  very  convenient  structure.  It  allows  a  physical  interpretation  of  the  data 
assimilation  process  in  terms  of  a  pseudo  scalar  potential  field  that  is  generated  by 
difference  between  current  measurements  and  predictions  of  those  measurements 
based  on  the  currently  estimated  vorticity  state.  This  pseudo  field  creates  an  additional 
velocity  field  that  causes  the  vorticity  state  to  adapt  to  the  measurements  at  the  same 
time  that  it  executes  the  normal  hydrodynamics. 

Estimates  of  the  computational  requirements  for  2D  and  3D  implementations  have 
been  made  and  are  described  in  Figure  3.1.3.  The  algorithm  is  relatively  efficient  and 
computational  loads  relatively  modest  for  restricted  grid  dimensions. 
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The  computation  using  the  Lagrangian  mode  has  been  implemented  in  a  2D  numerical 
simulation.  The  basic  equations  (in  3D  form)  used  are  summarized  in  this  section. 

The  computation  tracks  the  trajectories  of  a  cloud  of  discrete  vortex  elements  in  a  three 
dimensional  space.  The  position,  orientation,  and  strength  of  each  vortex  are  updated 
at  each  time  interval.  The  velocity  has  two  contributions: 

=  u(x)  +  VF(x,a)  (12) 

at 

The  first  (u)  is  the  physical  velocity  due  to  all  the  other  vortices  and  is  derived  as  the 
curl  of  a  local  vector  stream  function 

u  =  V  x  ¥  (13) 

where  the  stream  function  is  derived  by  solving  a  Poisson  equation  whose  source  is  the 
vector  vorticity 

V2vF  =  Q  (14) 

The  vorticity  is  computed  by  (vector)  summing  the  strengths  of  all  the  vortex  elements 
in  a  given  mesh  cell  and  dividing  by  the  mesh  cell  volume. 

The  second  velocity  term  is  the  gradient  of  a  scalar  potential,  evaluated  at  the  vortex 
location  which  arises  from  the  disparity  between  the  predicted  and  currently  observed 
data.  This  scalar  potential  is  the  dot  product  of  the  vortex  strength  and  an  error  vector 
potential  function  vector: 

F(x,a)  =  a-A,(x)  (15) 

The  error  vector  potential  As(x)  is  a  convolution  of  the  error  covariance  function  P(x) 
with  a  pseudo-stream  function  field  %  which  is  generated  from  the  measured  data: 

At{x)=  P(x)*x¥i(x)/RAXASAt  (16) 

In  this  equation  the  asterisk  represents  a  spatial  convolution.  The  error  covariance 
function  P(x)  arises  in  the  estimation  process  as  an  estimate  of  the  uncertainty  of  the 
estimated  locations  of  the  vortex  elements,  AXAt  and  AS  are  volume  elements  in  the 
state  and  measurement  space  respectively,  and  R  is  the  expected  measurement  error 
(variance)  of  the  line-of-sight  velocities. 

The  pseudo-stream  function  is  the  solution  of  a  vector  Poisson  equation  for  which  a 
pseudo  vorticity  field  Qs  is  the  source: 

V2vFs  =  Qs  (17) 

The  pseudo-vorticity  itself  is  calculated  from  the  measured  data  as  the  curl  of  a  vector 
field  consisting  just  of  the  observed  line-of-sight  velocities: 
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Q,(x)  =  VxZ,(x) 


(18) 


Here  the  vector  field  Zs(x)  is  constructed  from  the  lidar  measurements  of  the  line-of- 
sight  velocity  Zs(x,n)  as  a  sum  over  all  lines  of  sight  according  to 

z.(*)  =  XnZ*(x’n)  (19) 

where  n  is  a  unit  vector  in  the  direction  of  the  given  line  of  sight. 


Figure  3.1.  Kalman  Filter  Update  Sequence  :  Euierian 
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Calculate  net  velocity  at  each  vortex,  advect  all  vortices 


Calculate  true  velocity  at 
each  marker,  advect  all  markers 


Display  markers 


Calculate  new  true  vorticity  (x) 


t  =  t  +  dt 


Figure  3.2.  Kalman  Filter  Update  Sequence  :  Lagrangian 
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Update  Total  =  72  Ndata  +  8  N  Ig(N)  +  12  N 


Display  markers 


t  =  t  +  dt  : 


TOTAL  flops  (M  &  A)  per  cycle  :  102  Nv  +  72  Ndata  +  16  N  lg(N)  +  42  N  +  3 
Nv  =  500  vortices,  Ndata  =  16x64  =  1024  data  points,  N  =  64x16  =  1024  mesh  points  requires 
51  +  70  +  1 60  +  42  =  323  flops/meshpoint/cycle 


Figure  3.3  .  Arithmetic  operations  count  for  the  Lagrangian  mode  (2D  code) 
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These  equations  can  be  used  to  adapt  a  preexistent  field  of  vortex  elements  to  a  set  of 
measurements  while  maintaining  consistency  with  the  mass  and  momentum  equations 
of  motion  of  the  atmospheric  fluid. 

In  addition  to  moving  the  vortex  elements  to  better  fit  the  observations,  the  Kalman  filter 
can  also  change  their  strengths  and  orientations.  The  growth  rate  of  the  circulation  of 
vortex  element  p  is  proportional  to  the  negative  of  the  pseudo-stream  function 


da, 

dt  RAS&t  (20) 

where  aa  is  the  expected  rms  uncertainty  of  the  estimate  of  the  vortex  strength  ap.  A 
detailed  derivation  of  these  relationships  is  provided  in  the  Appendix. 


3.2  WIND  VELOCITY  ESTIMATION:  SIMULATIONS  AND  EXPECTED 
PERFORMANCE 

A  previously  developed  2D  simulation  hydro  code  has  been  modified  to  use  these 
Lagrangian  mode  equations.  Here  the  vorticities  and  stream  functions  have  only  a 
single  component  in  the  direction  perpendicular  to  the  simulation  plane. 

The  simulation  is  used  in  two  steps.  First  a  calculation  with  no  input  data  is  carried  out 
starting  from  some  selected  initial  vorticity  distribution.  The  radial  velocity  that  would 
be  sensed  along  selected  lines  of  sight  is  recorded.  This  recorded  data,  modified  by 
the  lidar  response  function  and  with  added  noise,  is  then  used  as  data  input  for  a  data 
assimilation  run. 

For  a  demonstration  of  the  process  we  have  constructed  the  velocity  field  of  a 
simulated  2D  wind  gust.  Here  the  initial  condition  consists  of  a  downward  moving 
parcel  of  air  initially  centered  at  an  a  altitude  of  about  700  meters  (Figure  1.2). 
Because  of  the  periodic  boundary  conditions  imposed  in  the  horizontal  direction,  there 
are  actually  a  row  of  wind  gusts  each  separated  by  the  width  of  the  computed  field. 
These  gusts  impinge  on  the  ground  surface  at  an  angle  of  about  45  degrees  and  each 
split  into  two  parts  moving  in  opposite  directions.  These  parcels  travel  horizontally  until 
they  encounter  other  parcels  entering  the  computational  mesh  from  adjacent  periodic 
image  cells.  The  collision  of  these  parcels  then  results  in  an  upward  travelling  jet. 

Simulated  lidar  data  is  recorded  along  each  of  the  lidar  lines  of  sight  shown  in  Figures 

1.2  and  1.3.  The  data  is  derived  by  convolving  the  calculated  line  of  sight  component  of 
the  fluid  velocity  with  the  lidar  range  response  function,  and  adding  simulated  photon 
and  speckle  noise.  This  data  is  recorded  and  is  to  be  input  to  a  second,  data 
assimilation  code. 
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Implementation  of  the  Kalman  Filter 


Computational  sequences  for  Eulerian  and  Lagrangian  modes  are  outlined  in  Figures 
3.1  and  3.2.  In  these  figures  the  greyed  areas  represents  the  update  that  results  from 
the  incorporation  of  the  measured  data,  the  remainder  represents  the  calculation  of  the 
natural  inviscid  hydrodynamic  motion.  In  both  modes  the  calculation  cycle  starts  with 
the  three  dimensional  distribution  of  vector  vorticity.  At  the  initial  time  this  distribution 
can  be  selected  arbitrarily,  and  is  typically  set  to  yield  a  motionless  or  uniformly  moving 
atmosphere.  This  vorticity  serves  as  the  source  function  for  a  Poisson  equation  given 
in  Eq.  (14).  Gradients  (the  curl)  of  the  stream  function  solution  are  evaluated  at  the 
location  of  individual  vortex  threads  to  yield  the  fluid  advection  velocity.  In  the  absence 
of  measurement  input  this  velocity  is  used  to  advect  the  vortex  threads  with  the  local 
fluid  for  the  given  time  interval.  A  revised  vorticity  distribution  is  then  calculated  based 
on  the  new  vortex  positions  and  the  process  repeated. 

When  measurements  from  the  lidar  are  to  be  incorporated,  the  line-of-sight  radial 
velocity  is  evaluated  at  each  of  the  range  gates  using  the  current  velocity  field 
predictions.  The  difference  between  these  predictions  and  the  measured  values  are 
inserted  into  an  x-y  grid  at  the  measurement  locations  and  the  rest  of  the  grid  is  set  to 
zero.  The  Eulerian  and  the  Lagrangian  procedures  differ  primarily  in  the  way  these 
differences  are  then  assimilated. 

In  the  Eulerian  calculation  mode  the  3D  data  error  field  is  passed  through  a  spatial  filter 
to  yield  an  increment  of  vorticity  at  each  point  in  the  spatial  grid.  This  vorticity 
increment  field  is  then  added  to  the  current  vorticity  estimate  produced  by  the  most 
recent  hydrodynamic  prediction  step  of  the  calculation.  A  new  set  of  discrete  vortex 
elements  is  then  built  to  represent  this  updated  vorticity  field  and  the  process  repeated. 

In  the  Lagrangian  mode  the  data  errors  for  different  look  directions  are  combined  to 
form  a  3D  vector  error  field  which  consists  just  of  the  vector  difference  between  the 
measured  line  of  sight  components  of  the  wind  field  and  their  predictions.  Only 
locations  and  components  that  actually  contribute  to  the  measurements  contribute  to 
this  vector  field;  all  other  velocity  values  are  set  equal  to  zero.  The  curl  of  this  psuedo- 
wind  field  is  then  used  as  the  source  function  for  a  Poisson  equation. 

The  solution  of  this  Poisson  equation,  after  being  passed  through  a  second  spatial 
filter,  is  used  for  two  purposes:  First,  it  serves  a  source  for  generating  new  vorticity.  To 
this  point  the  Lagrangian  calculation  is,  in  principle,  equivalent  to  the  update  process  in 
the  Eulerian  mode.  In  the  Lagrangian  mode,  however,  the  Poisson  solution  field  is  also 
used  to  create  a  scalar  pseudo-potential  whose  gradients  are  added  to  the  physical 
velocity  field  and  the  total  used  to  advect  the  discrete  vortex  elements.  This  allows  two 
methods  for  the  correction:  one  is  the  creation  of  new  vorticity,  the  other  is  the 
advection  of  existent  vorticity. 

Since  the  pseudo-potential  is  equal  to  the  projection  of  the  vortex  strength  on  the 
vector  solution  of  the  Poisson  equation,  vortex  elements  having  opposite  strengths  will 
move  in  opposite  directions  in  this  error  correction  field.  This  latter  advection  allows  a 
preexisting  vorticity  field  to  adapt  to  measurement  input  without  necessarily  requiring 
the  creation  of  new  vorticity. 
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Figure  1.2  shows,  as  a  function  of  time,  the  developing  inferred  vector  field  produced 
by  the  Kalman  filter  for  a  wind  gust  simulation.  Figure  1.2a  shows  the  velocity  state 
shortly  after  initiation  of  the  measurement.  In  this  example  the  wind  field  was  created 
by  a  simple  vortex  pair  representation  of  a  wind  gust.  Data  has  been  acquired  along 
three  lines  of  sight.  In  the  reconstruction,  the  inferred  atmosphere  is  assumed  to  be 
initially  motionless.  The  Kalman  filter  then  generates  and  advects  vorticity  to  allow  the 
atmosphere  to  adapt  to  the  measurements. 

The  upper  panels  of  each  of  Figures  1.2a  and  1.2b  show  the  inferred  wind  vectors  and 
the  locations  of  the  lines  of  sight.  The  six  lower  panels  in  each  of  Figures  1.2a  and  b 
show  the  simulated  measured  data  (solid  lines)  and  the  current  predictions  (dots). 
Panels  containing  asterisks  are  the  measured  lines  of  sight.  The  intermediate  panels 
show  true  and  predicted  data  along  intermediate  sight  lines.  At  the  time  of  this  first 
display  there  are  significant  differences  between  the  predictions  of  the  radial  winds 
and  the  measurements  along  both  participating  and  non-participating  lines  of  sight. 
The  difference  between  these  values  drives  the  adaption  process.  The  adaption  that 
has  taken  place  at  this  time  is  due  primarily  to  the  imposition  of  the  continuity 
constraint.  With  time,  both  the  momentum  constraint  and  the  continuity  constraint 
continue  to  cause  the  solution  field  to  adapt  to  the  data.  Later  in  time  (Figure  1 .2b) 
these  differences  are  reduced  to  the  point  that  the  predictions  are  effectively  identical 
to  the  measurements.  At  this  time  the  modeled  atmospheric  field  is  fully  consistent  with 
the  data  as  well  as  with  all  the  hydrodynamic  constraints  embodied  in  the  fluid 
equations  of  motion.  At  this  later  time  the  predicted  data  field  also  agrees  very  well 
with  true  radial  velocity  on  the  non-participating  sight  lines  indicating  an  accurate 
reconstruction  of  the  entire  vector  field. 

Reconstructions  of  a  more  complex  wind  field  are  shown  in  Figure  1.3.  The  true 
velocity  field  at  time  101  seconds  is  shown  in  panel  a  and  reconstructions  using  3  and 
6  lines  of  sight  in  panels  b  and  c.  In  this  case  the  reconstruction  with  3  lines  of  sight  is 
less  satisfactory  even  though  the  deduced  field  at  the  measurement  points  is  very 
accurate.  With  6  lines  of  sight  the  inferred  field  is  well  predicted  throughout  the 
measurement  field. 

Examination  of  these  simulations  can  be  used  to  establish  the  ability  of  the  retrieval 
process  to  reconstruct  the  wind  field  solely  from  radial  wind  data  as  the  scan  strategy  is 
changed.  In  general  the  process  appears  to  converge  to  reasonable  approximations  of 
the  true  flow  field  if  sufficient  number  of  lines  of  sight  are  chosen  in  the  vertical  plane. 

In  the  Kalman  filter  analysis  as  presently  implemented,  both  the  constraints  of 
incompressibility  and  momentum  conservation  aid  the  retrieval  of  the  vector  wind  field 
from  the  scalar  radial  wind  measurements.  The  momentum  conservation  constraint  is 
basically  a  method  for  tracking  the  advection  of  persistent  patterns  of  radial  velocity. 
Although  not  shown  in  these  simulations,  detection  of  patterns  of  aerosol  scattering 
intensity  can  also  be  used  to  assist  the  velocity  retrieval. 

A  similar  ability  to  detect  and  track  patterns  of  absorption  of  water  vapor  can  be 
included  by  incorporating  absorption  into  the  measurement  model  and  allowing  for 
passive  advection  of  water  vapor  in  the  state  description. 
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Computation  Issues 


Figure  3.3  shows  an  approximate  count  of  floating  point  operations  for  the  Lagrangian 
mode  as  implemented  in  the  2D  simulation.  One  multiply  and  add  is  counted  as  one 
floating  point  operation.  These  numbers  were  evaluated  by  counting  the  operations  in 
the  coded  2D  calculation  used  to  produce  the  simulations  shown  in  this  report.  The 
count  excluded  all  graphics  output  as  well  as  any  I/O  operations.  For  the  sparse  lidar 
sampling  selected  (6  lines  of  sight)  and  the  modest  number  of  vortices  allowed  (on  the 
average  a  given  mesh  cell  is  populated  50%  of  the  time)  the  computation  is  dominated 
by  the  Fourier  transforms  required  to  generate  the  real  and  pseudo  velocity  fields.  For 
this  coding  somewhat  more  than  300  operations  per  spatial  mesh  point  per  cycle  are 
required.  For  the  64x16  mesh,  this  corresponds  to  300,000  operations  per  cycle 
expected  for  the  2D  simulation.  In  3D  the  number  of  operations  per  mesh  point  is 
roughly  tripled,  since  three  components  of  vorticity  and  streamfunction  are  required 
rather  than  just  the  one  needed  in  2D.  Thus,  for  a  64x32x16  mesh,  we  anticipate  a 
computation  load  of  order  30  million  multiply  and  add  operations  per  cycle.  Given  that 
cycle  times  of  one  to  several  seconds  are  anticipated,  the  present  SkyBolt  board  used 
in  CTI's  signal  processor  should  be  (just  barely)  capable  of  handling  this  load  in  a  real 
time  operational  mode. 


4.0  PRELIMINARY  VELOCITY  MEASUREMENTS 

The  2.1  micron  coherent  lidar  system  is  currently  being  integrated  into  a  mobile  lidar 
van  (see  Section  6.0).  Prior  to  its  installation,  the  system  has  operated  out  of  a  building 
at  CTI's  Table  Mountain  test  facilities.  The  system  has  a  10  cm  aperture  and  uses  a 
computer-controlled,  super-hemispherical  scanner  to  produce  Plan  Position  Indicator 
(PPI)  and  Range  Height  Indicator  (RHI)  displays.  The  system  transmits  20  mJ/pulse 
with  a  pulse  repetition  frequency  (PRF)  of  5  Hz  and  a  full  width  half  maximum  (FWHM) 
pulsewidth  of  roughly  200  nsec  (range  resolution  -30  meters).  The  200  nsec 
pulsewidth  permits  the  type  of  high  resolution  velocity  measurements  required  to 
understand  the  dynamics  of  the  turbulent  boundary  layer. 

Signal  returns  are  digitized  at  a  rate  of  100  megasamples/sec  and  processed  with  CTI’s 
Real  Time  Lidar  Processor  (RTLP)  and  displayed  in  a  variety  of  formats.  Figures  1.7a 
and  1.7b  show  two  examples  of  wind  velocity  measurements  made  with  the  system. 

Figure  1.7a  is  a  plot  of  the  signal  frequency  content  (spectrum)  versus  range  along  a 
single  line  of  sight.  The  lidar  was  pointed  slightly  above  horizontal  (+3°  elevation  angle) 
and  to  the  north  of  Table  Mountain  toward  a  hillside.  In  each  of  64  equally  spaced 
range  gates,  a  64  point  Fast  Fourier  Transform  (FFT)  is  computed.  The  sample 
spacing  is  10  nsec  (1.5  m)  such  that  the  width  of  each  range  gate  is  96  m.  The  FFT 
amplitudes  are  color-coded  and  displayed  as  a  horizontal  'stripe'  along  the  frequency 
axis  in  Figure  1.7a.  This  figure  represents  one  second  of  data  (5  pulses).  The  range 
gates  (y-axis)  cover  the  line  of  sight  range  window  from  1  km  to  5.54  km.  The 
frequency  axis  (x-axis)  extends  from  -15  MHz  to  +15  MHz  (or  +15  m/sec  to  -15  m/sec  at 
the  2  micron  wavelength).  A  frequency  of  0  MHz  corresponds  to  no  Doppler  shift  and 
thus  implies  a  zero  radial  (relative  to  the  lidar)  velocity.  The  color  lookup  table  at  the 
bottom  of  the  screen  extends  from  10  (blue)  to  70  (red)  in  arbitrary  log  scale  units  (a 
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total  range  of  60  dB).  The  red  regions  indicate  the  location  in  frequency  and  range  of 
the  energetic  portions  of  the  signal  returns.  The  blue  regions  indicate  the  absence  of 
signal  and  are  due  to  the  local  oscillator  shot  noise  floor.  At  a  range  near  5  km,  an 
extremely  strong  return  centered  about  0  MHz  is  due  to  the  distant  hillside.  The  0  MHz 
frequency  is  expected  as  the  hillside  is  not  moving.  Beyond  the  hillside,  no  signal 
radiation  is  available  for  detection,  and  the  blue-colored  noise  floor  is  all  that  remains. 
In  front  of  the  hillside,  the  aerosol  returns  permit  accurate  wind  velocity  estimation.  A 
majority  of  the  aerosol  returns  are  centered  near  -6  MHz  and  indicate  a  radial  wind 
velocity  of  +6  m/sec  (positive  velocity  away  from  the  lidar).  Note  the  interesting  velocity 
structure  within  the  first  2.5  km,  especially  the  5  MHz  (5  m/sec)  shear  from  1.4  km  to 
1.9  km. 

In  Figure  1.7b,  we  plot  the  FFT  first  moment  (velocity)  versus  range  in  a  rectangular, 
side-looking  PPI  display.  For  this  plot,  the  lidar  scanned  through  90°  in  azimuth  with 
the  elevation  angle  held  constant  at  +3°.  Please  note  that  this  data  was  collected  on  a 
different  day  from  that  in  Figure  1.7a  The  y-axis  is  the  azimuth  angle  and  extends  from 
70°  to  160°  (90°  is  due  North)  and  the  x-axis  is  range  and  extends  from  1  km  to  5.2  km. 
Radial  velocities  are  computed  from  the  amplitude  thresholded  first  moment  of  64-point 
FFTs  at  each  of  64  equally-spaced  range  gates.  The  FFTs  from  five  adjacent  lines  of 
sight  are  combined  to  compute  each  mean  velocity  output  record  (five  pulse 
averaging).  Each  output  record  is  plotted  as  a  horizontal  stripe  of  radial  wind  velocity 
versus  range.  The  y-axis  location  for  each  horizontal  stripe  corresponds  to  the 
averaged  azimuth  angle  for  the  five  lines  of  sight.  The  color  lookup  table  at  the  bottom 
of  the  screen  extends  from  -8  m/sec  (blue)  to  +8  m/sec  (red).  Positive  velocities  are 
again  defined  as  those  away  from  the  lidar.  At  an  azimuth  angle  of  150,  a  near  Field 
telephone  pole  blocks  the  beam,  eliminates  aerosol  returns,  and  results  in  poor  velocity 
estimates.  Over  the  90°  azimuth  scan,  the  radial  wind  velocity  shifts  from  values  near 
+3  m/sec  to  values  near  -7  m/sec.  From  an  azimuth  angle  of  90°  to  an  azimuth  angle 
of  135°,  interesting  small-scale  (1  to  2  m/sec  over  ranges  of  200  to  400  m)  wind  eddies 
are  clearly  present. 


System  Performance  and  SNR  Model  Validation 

Performance  model  validation  requires  a  large  data  base  of  well-documented  field 
measurements.  Two  of  the  more  important  model  inputs/assumptions  are  the  system 
SNR  and  the  degree  (accuracy)  to  which  the  vector  velocities  can  be  inferred  from 
radial  velocity  measurements.  The  degree  to  which  the  performance  prediction  models 
can  be  validated  greatly  impacts  the  prototype  system  design. 

The  validity  of  the  lidar  system  SNR  model  relies  heavily  on  how  well  atmospheric 
parameter  dependencies  at  the  operating  wavelength  are  known  both  horizontally  and 
vertically.  The  atmospheric  parameters  of  interest  are  the  volume  backscatter 
coefficient,  the  atmospheric  extinction  coefficient,  and  the  refractive  turbulence 
structure  coefficient.  In  addition,  the  system  efficiency  (photodetector  quantum 
efficiency,  optical  efficiency,  and  heterodyne  mixing  efficiency)  must  be  known  fairly 
accurately. 
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Figure  4. 1  Receiver  power  vs  range 
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As  a  first  step  toward  system  SNR  model  validation,  under  a  contract  to  Wright 
Laboratories,  CTI  recently  completed  field  measurements  of  SNR  versus  range  and 
elevation  angle.  The  2.09  mm,  21  mJ,  5  Hz,  10  cm  aperture  (8  cm  e'2  intensity 
diameter)  system  was  employed.  The  measured  values  of  SNR  were  compared  to 
those  predicted  by  the  model.  Table  1  shows  an  example  of  the  results.  The  skies 
were  overcast  during  the  measurements,  so  that  refractive  turbulence  effects  should  be 
small  (although  no  direct,  quantitative  measurement  of  refractive  turbulence  was 
made).  For  comparison  purposes,  model  predictions  for  the  standard  refractive 
turbulence  structure  coefficient  model  (Hufnagel)  are  given  in  addition  to  this  same 
model  reduced  by  a  factor  of  0.01.  As  expected,  the  latter,  'weak'  turbulence  model 
predictions  better  match  the  measured  values.  The  agreement  is  quite  good  for  all 
elevation  angles  and  all  ranges. 

Because  the  lidar  system  efficiency  is  currently  not  calibrated,  Table  1  says  nothing 
about  the  absolute  value  of  the  volume  backscatter  coefficient.  Rather,  it  shows  that 
the  product  of  system  efficiency  and  atmospheric  backscatter  (hb)  matches  well  to  that 
employed  by  the  model  (0.1  x  1x10*®  m^sH  at  ground  level).  More  importantly, 
however,  Table  1  shows  that  the  parametric  system  dependencies  (i.e.,  SNR  versus 
altitude)  are  well-modeled  by  the  simulation. 

Receiver  power  versus  range  data  taken  between  December  and  June  1992  has  been 
used  to  produce  Figure  4.1.  This  figure  shows  the  time  history  of  the  maximum 
horizontal  and  vertical  wind  measurement  ranges.  We  define  the  maximum 
measurement  range  as  that  range  where  the  narrowband  SNR  drops  to  0  dB  (the 
average  signal  contribution  to  the  detector  output  is  equal  to  the  average  shot  noise 
level).  Through  pulse  averaging  of  three  to  four  pulses  to  reduce  noise  speckle 
fluctuations,  a  narrowband  SNR  of  0  dB  or  more  typically  results  in  accurate  velocity 
measurements. 

For  the  horizontal  data  of  Figure  4.1,  the  lidar  beam  was  aimed  north  at  an  elevation 
angle  of  -0.4°  which  allows  the  beam  to  pass  just  above  the  mountains  to  the  north  of 
Table  Mountain.  For  the  vertical  data  of  Figure  4.1,  the  beam  was  aimed  approximately 
straight  up  (elevation  angle  of  90°).  For  all  the  data,  the  transmitted  beam  was 
collimated  and  the  transmitted  pulse  energy  was  roughly  20  mJ.  The  horizontal  data 
comprises  73  irregularly  spaced  measurements  (due,  for  example,  to  system  shutdown 
for  modifications  and  upgrades).  The  maximum  horizontal  measurement  range 
fluctuates  from  6  km  to  30  km.  Sixty-eight  percent  of  the  data  is  >  15  km  and  87%  is  > 
10  km.  The  vertical  data  represents  the  maximum  measurement  range  obtained  in  the 
absence  of  clouds.  No  vertical  data  was  taken  when  low-level  clouds  interfered  with 
the  measurement.  The  vertical  data  comprises  40  measurements  that  vary  from  1.3  km 
to  9  km.  Seventy-five  percent  of  the  data  is  >  3  km  and  87%  is  >  2  km. 

Since  correlations  in  aerosol  backscatter,  atmospheric  refractive  turbulence, 
atmospheric  transmission,  and  lidar  system  performance  have  not  been  made,  the  data 
is  only  intended  to  show  the  time  history  of  the  maximum  measurement  range  for  the 
ground-based  2  mm  system  in  Boulder,  CO.  Variations  in  atmospheric  refractive 
turbulence  and  aerosol  backscatter  can  cause  large  fluctuations  in  SNR  and  probably 
large  contributors  to  the  observed  measurement  range  fluctuations. 
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5.0  LASER  TECHNOLOGY  CONSIDERATIONS 


5.1  BACKGROUND 

CTI  has  led  the  development  of  efficient  diode-laser  pumped  2  micron  lasers  based 
around  TM-doped  YAG.  Figure  5.1  shows  data  from  a  relatively  early  experiment. 
Here  pulse  energy  is  plotted  against  pulse  spacing,  i.e.,  the  inverse  PRF.  Using  only  6 
W  of  pump  power,  this  diode-pumped  laser  produced  4.5  mJ  at  100  Hz,  and  3  mJ  at 
200  Hz.  Later  improvements  yielded  >  5  mJ  near  100  Hz  and  4.2  mJ  at  200  Hz.  As 
shown  in  Section  2  these  characteristics  are  close  to  those  needed  for  a  coherent  lidar 
turbulence  probe  for  the  atmospheric  boundary  layer. 

One  of  the  prime  areas  of  commercial  interest  for  lidar  systems  is  for  airborne 
windshear  detection.  Based  upon  the  laser  results  discussed  above,  CTI  is  presently 
collaborating  with  a  major  aerospace  company  to  develop  a  flightworthy  windshear 
detector.  At  present  this  transceiver  is  scheduled  for  delivery  in  April  1993. 

Figure  5.2  shows  a  sketch  of  the  transceiver  as  it  is  being  developed  .  It  consists  of 
two  lasers,  the  master  oscillator  (MO)  and  the  slave  oscillator  SO).  The  MO  power  is 
split  into  several  signals  using  optical  fibers/couplers.  Part  of  the  MO  signal  is  injected 
into  the  SO  which  forces  this  Q-switched  laser  to  oscillate  in  a  single  longitudinal  mode. 
The  SO  output  is  sent  through  the  transmit/receive  (T/R)  switch  out  through  a  telescope 
and  scanner.  Radiation  is  collected  with  the  same  aperture,  and  is  coupled  into  an 
optical  fiber,  and  mixed  with  another  port  of  the  MO  power.  That  signal  is  detected  by  a 
photodetector  (PD1),  which  is  connected  to  the  signal  processor.  To  ensure  minimum 
uncertainty  in  the  velocity  estimates,  a  second  detector  (PD2)  is  used  to  monitor,  on  a 
pulse-to-pulse  basis,  the  exact  frequency  shift  between  MO  and  SO. 

In  parallel  with  the  hardware  development,  CTI  is  also  working  towards  increasing  laser 
pulse  energies.  We  are  presently  able  to  achieve  very  high  brightness  from  optical 
fibers  pumped  with  the  output  from  several  separate  diodes.  Fiber  pumping  allows  one 
to  separate  the  pump  optics  from  the  laser  rod.  Using  three  diodes  coupled  into  a 
single  200  micron  diameter  fiber,  and  a  fourth  diode  coupling  directly  into  the  laser  rod, 
we  have  achieved  CW  output  powers  exceeding  2  W.  Much  higher  values  are 
expected  shortly,  and  Q-switched  energy  increases  of  almost  a  factor  of  two  are  also 
expected.  For  further  increases  in  energy,  more  pump  diodes  can  relatively  easily  be 
added.  It  is  also  important  to  note  that  2  micron  lasers,  such  as  Tm:YAG,  are  tunable 
over  quite  large  ranges.  The  CTI  MO  laser  has  been  continuously  tuned  from  about 
2007  nm  to  2023  nm  (corresponding  to  >  1  THz),  and  the  SO  has  been  injection- 
seeded  to  2021  nm.  Since  the  emission  peak  occurs  at  about  2012  nm,  this  represents 
a  detuning  of  9  nm.  The  great  tuning  range,  presently  limited  by  the  tuning  etalons  in 
the  MO,  is  important  since  it  enables  one  to  select  an  optimal  operating  wavelength 
quite  freely. 


5.2  LASER  TECHNOLOGY  FOR  DIAL  MEASUREMENTS 

DIAL  meaurements  of  the  type  discussed  above  require  a  high  prf  source  (typically  300 
Hz)  that  can  be  switched  between  an  absorbing  (A)  and  non-absorbing  (N)  wavelength. 


44 


45 


Block  Diagram  of  Coherent  Laser  Radar  System 

Coherent  Technologies,  Inc. 
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Figure  5.2.  Windshear  Transceiver  Unit  Diagram 


Wavelength  switching  can  be  achieved  in  two  basic  ways.  Either  one  source  is  used  to 
tune  between  the  A  an  N  wavelengths,  or  two  lasers  are  used,  and  the  transmitter  is 
switched  between  lasers.  The  first  option  is  the  more  elegant,  but  also  has  some 
potential  porblems  that  require  fundamental  studies.  The  main  difficulty  is  not 
necessarily  achieving  rapid  tuning.  Rather  it  relates  to  ensuring  that  the  laser 
wavelength  settles  down  to  the  desired  ones  quickly,  that  repeatability  is  high,  and  that 
the  laser  always  remains  single-frequency.  The  second  option  is  more  straightforward, 
since  two  separate  lasers  can  be  stabilized  quite  well  in  wavelength. 

We  propose  that  the  simplest  and  most  cost-effective  system  is  a  hybrid  one.  Two  CW 
master  oscillators,  stabilized  to  different  wavelengths,  are  used  to  injection-seed  a 
single  slave  oscillator.  The  main  advantage  of  this  approach  is  that  only  one  slave 
oscillator  is  needed,  which  saves  considerable  cost  and  ensures  that  the  energies 
transmitted  at  the  two  wavelengths  are  the  same.  The  output  energy  of  an  injection- 
seeded  laser  is  independent  of  the  seed  power,  so  small  power  differences  between 
the  two  master  oscillators  are  not  important. 

Wavelength  Choices 

In  selecting  suitable  wavelengths  for  DIAL  measurements,  several  issues  must  be 
considered.  First,  the  non-absorbing  one  should  have  quite  a  low  attenuation  over  the 
desired  range.  Second,  the  absorbing  one  should  have  an  attenuation  that  is  not  too 
weak,  but  also  not  too  strong.  If  it  is  weak,  the  differential  absorption  becomes  more 
difficult,  while  if  it  is  too  strong,  there  may  not  be  enough  signal  from  the  farthest 
desired  range.  The  A  wavelength  must  also  be  directly  associated  with  water 
absorption,  and  should  ideally  not  be  contaminated  by  absorption  due  to  other  species. 
Figure  1.4a  shows  the  MLS  attenuation  spectrum  in  the  region  of  tuning  of  Tm:YAG 
lasers.  There  is  a  clear  C02  overtone  band  covering  much  of  the  region,  as  well  as 
other  peaks  associated  with  water  and  N02.  Figure  1.4b  shows  the  water  spectrum 
alone.  Several  of  the  water  lines  overlap  directly  with  C02  lines  and  are  therefore 
unsuitable  for  this  application.  There  are  also  several  lines  that  do  not  overlap,  notably 
the  one  marked  with  an  arrow  and  the  letter  A  at  2016.974  nm  (vacuum  wavelength  - 
the  air  wavelength  will  be  approximately  0.5  nm  shorter).  Detailed  examination  of  that 
line  shows  an  extinction  coefficient  of  3.389  km  (29  dB/km  roundtrip),  and  that  the  full 
width  at  the  95%  peak  value  is  approximately  0.019  nm,  corresponding  to  1.41  Hz.  The 
width  is  important,  as  it  directly  relates  to  how  frequency-stable  the  laser  source  must 
be.  Further  examination  of  the  peak  shows  that  neighboring  C02  lines  contribute  <  3% 
of  the  total  attenuation.  This  line  is  therefore  an  excellent  choice  for  short-range  water 
vapor  measurements. 

For  the  non-absorbing  wavelength,  the  best  choices  is  to  operate  at  2017.713  nm.  The 
valley  shows  an  adsorption  of  0.1  per  km  at  that  point  (0.4  db  roundtrip  attenuation  per 
km).  The  width  of  the  valley  floor  is  approximately  0.078  nm,  corresponding  to  a 
frequency  width  of  5.8  GHz. 

Both  of  these  wavelengths  are  easily  accessible  with  an  injection-seeded  Tm:YAG 
laser.  Continous  tuning  from  2006  to  2023  nm  has  been  demonstrated. 
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Wavelength  Stabilization 


Stabilization  of  the  two  master  oscillators  (referred  to  as  MO-N  and  MO-A)  c^n  oe  done 
relatively  simply.  In  CTI's  master  oscillator  design,  two  intracavity  etalons  are  used  to 
control  the  emission  wavelength.  Tuning  is  accomplished  by  rotating  the  etalons  or 
varying  their  temperatures.  The  thin  etalon  determines  the  emission  wavelength  .  while 
the  thick  one  ensures  that  the  laser  operates  in  a  single  transverse  mode.  The  mode 
spacing  is  determined  by  the  free  spectral  range  of  the  laser  cavity,  approximately  1.5 
GHz  in  our  case.  The  thick  etalon  is  the  element  most  sensitive  to  temperature 
variations,  since  it  temperature  tunes  at  a  rate  of  approximately  1.4  GHz  per  degree  C. 
Our  lasers  have  the  thick  etalon  stabilized  to  0.1  degrees  C,  which  ensures  that  the 
wavelength  drift  is  no  more  than  about  140  MHz.  This  has  experimentally  been 
determined  to  be  sufficient  to  keep  the  laser  from  oscillating  on  more  than  one 
longitudinal  mode.  Note  that  preventing  simultaneous  oscillation  at  several  frequencies 
is  the  real  reason  for  keeping  the  laser  temperature  stabilized  .  Atmospheric  features 
are  typically  broad  compared  with  140  MHz,  and  hence  absolute  wavelength  accuracy 
to  that  level  is  not  required. 

Given  these  parameters,  operation  at  a  particular  wavelength  is  really  a  question  of 
settability.  Once  the  oscillating  wavelength  has  been  set,  simple  temperature  control  of 
the  laser  will  keep  it  at  the  right  wavelength.  Note  here  that  the  absolute  wavelength 
stability  requirement  is  1.4  GHz  for  the  A  wavelength,  and  almost  6  GHz  for  the  N 
wavelength.  The  only  real  difficulty  is  in  setting  the  wavelengths  to  their  proper  values 
from  the  outset.  A  good  monochromator  can  be  used  to  set  the  wavelength  to  will 
below  1  nm,  which  is  about  one  order  of  magnitude  away  from  the  desired  accuracy. 
One  way  of  getting  to  the  0.019  nm  level  is  to  use  a  wavemeter.  A  second  method  is 
simply  to  tune  the  laser(s)  and  look  for  maximum  (N)  and  minimum  (A)  signals.  This  is 
probably  a  more  accurate  method,  since  it  directly  probes  the  absorption  features  we 
are  interested  in  measuring. 

While  we  do  not  anticipate  difficulties  in  setting  the  laser  wavelengths  properly,  we  also 
note  that  there  are  active  methods  available  for  ensuring  frequency  stability.  To  do  this 
one  can  build  a  water  vapor  absorption  cell,  through  which  the  MO-A  beam  is  passed. 
With  a  sufficient  water  vapor  concentration  and  a  long  enough  path,  the  transmitted 
power  will  vary  dramatically  as  the  laser  wavelength  is  tuned  through  the  absorption 
resonances.  We  estimate  that  aim  long  cell  maintained  at  40  degrees  C,  and 
saturated  with  water  vapor,  through  which  the  beam  makes  a  single  pass,  will  show  a 
12  %  change  in  absorption  at  the  peak,  compared  to  off  the  peak.  Such  a  cell  is 
straightforward  to  construct,  and  furthermore,  by  letting  the  beam  make  multiple 
passes,  or  increasing  the  cell  temperature,  will  increase  the  attenuation  more. 

The  idea  is  then  to  measure  the  cell  transmission  and  use  a  servo  to  lock  onto  the 
water  resonance.  Since  drit.s  of  the  MO  wavelength  are  expected  to  be  small  over 
several  second  timeframes,  one  of  the  MO  etalons  can  simply  be  temperature  tuned  to 
close  the  servo  loop.  It  is  also  possible  with  this  active  stabilization  scheme  to  move  off 
the  absorbing  peak  in  a  controllable  manner,  so  that  the  measurement  range  can  be 
increased. 

Switching  between  the  two  MOs  can  be  achieved  most  simply  by  using  a  mechanical 
chopper,  which  can  also  be  used  to  trigger  the  slave  oscillator  Q-switch.  For  example, 
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Figure  5.3  Frequency-Switched  DIAL  Transmitter 


a  10-slot  chopper  wheel  running  at  15  rps,  together  with  a  prism/beamsplitter  can  be 
used  to  gate  either  of  the  two  MO  beams.  With  a  simple  LED/detector  pair  mounted  to 
the  chopper,  low  and  high  transitions  can  be  used  to  trigger  the  Q-switch  to  fire.  The 
optical  assembly  also  has  two  output  ports,  so  that  one  port  can  be  used  to  seed  the 
slave  oscillator,  while  the  other  one  is  used  to  provide  a  local  oscillator  signal.  Figure 
5.3  shows  a  sketch  of  how  such  a  dual-MO  laser  would  be  configured. 


6.0  PROTOTYPE  MEASUREMENT  SYSTEM:  CONCEPTUAL  DESCRIPTION 

A  program  to  develop  a  lidar  boundary  layer  probe  as  described  in  this  report  would 
have  the  following  ingredients: 

-  construction  or  acquisition  of  a  1  or  2  watt  mean  power,  pulsed  diode  pumped 
laser  that  can  operate  at  several  hundred  Hz. 

-  incorporation  of  this  lidar  into  a  mobile  scanning  system 

-  development  or  acquisition  of  the  3D  hydrodynamic  model 

-  extension  of  the  Kalman  filter  software  to  3D 

-  verification  of  the  vector  wind  retrieval  process  using  LES  output  to  simulate 
radial  wind  data  for  realistic  boundary  layer  wind  fields 

-  construction  of  the  DIAL  multiple  wavelength  switching  capability 

-  acquisition  of  3D  radial  wind  data  at  a  suitably  instrumented  location  (airport) 

-  demonstration  of  vector  wind  retrieval 

-  acquistion  of  simultaneous  3D  wind  and  DIAL  data  at  an  instrumented  location 

-  demonstration  and  verification  of  water  vapor  retrieval 


Under  a  Wright  Laboratory  contract,  CTI  is  currently  installing  a  2  micron  coherent  laser 
radar  transceiver  into  a  compact,  mobile  van.  The  lidar  van  is  shown  in  Figure  6.1. 
The  lidar  van  includes  a  computer-controlled,  super-hemispherical  scanner  that  is  fully 
integrated  into  the  signal  processing  system.  The  lidar  van  is  temperature  controlled  to 
operate  in  all  environments. 

The  proposed  diode-pumped  2.0  micron  coherent  laser  radar  measurement  system  will 
be  smaller  than  the  existing  2.1  micron  system  and  could  be  easily  installed  on  the  lidar 
van  optical  bench.  Because  the  lidar  system  is  mobile,  boundary  layer  measurements 
could  be  made  in  a  variety  of  geographic  locations. 


50 


Real  Time  Lidar  Processor  (RTLP) 


The  Real  Time  Lidar  Processor  is  a  very  flexible,  completely  programmable,  data 
acquisition  and  analysis  workstation.  A  system  block  diagram  is  shown  in  Figure  6.2. 
The  system  uses  an  Analytek  waveform  sampler  for  data  acquisition,  SkyBolt  i860 
applications  accelerators  for  signal  processing,  a  SUN  SPARCstation  for  user  interface 
and  display,  and  other  standard  VME  modules  for  scanner  control  and  recording.  The 
Analytek  sample  module  uses  a  patented,  VLSI  analog  storage  module  for  burst 
acquisition  and  can  operate  at  programmable  rates  of  1  MHz  to  500  MHz,  in  1  MHz 
steps  with  12  bit  accuracy.  The  module  can  be  operated  in  either  single  (2  GHz 
maximum  sampling  rate)  or  dual  (1  GHz  maximum  sampling  rate)  channel  modes.  The 
high-speed  SkyBolt  DSP  is  an  application  accelerator  based  on  a  40  MHz  i860  capable 
of  80  MFLOPS  of  peak  performance,  and  an  i960  processor  to  manage  memory  and 
I/O.  The  SkyBolt  can  be  programmed  in  either  FORTRAN  or  C,  and  is  supported  with  a 
vast  library  of  vector  calls  capable  of  performing  at  or  better  than  50%  of  the  peak 
theoretical  performance. 

The  RTLP  performs  real-time  A-scope  and  range-resolved  periodogram  (spectral) 
processing  of  aerosol  data.  The  fully  integrated  system  provides  processing,  display, 
recording,  and  playback  capabilities  with  real-time,  color-coded  graphics.  The  RTLP  is 
fully  programmable  with  standard  outputs  including  plan  position  indicator  (PPI) 
displays,  color-coded  range  height  indicator  (RHI)  displays,  radial  velocity  mean  and 
width  versus  range  plots,  and  velocity  versus  height  plots  resulting  from  velocity 
azimuth  (VAD)  scans. 

The  RTLP  can  be  easily  upgraded  to  incorporate  the  requisite  software  for  combined  3- 
D  wind  velocity  and  water  vapor  measurement  computations.  Due  to  the  intense 
computational  loads  expected,  completely  real  time  operation  may  not  be  possible.  For 
example,  100  MFLOPS  of  sustained  performance  may  be  required.  In  such 
circumstances,  the  RTLP  would  collect  data  for  a  period  of  10  to  15  minutes,  process 
data  and  output  results  for  40  to  50  minutes,  process  more  data  for  10  to  15  minutes, 
process  data  and  output  results  for  40  to  50  minutes,  and  so  on. 
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APPENDIX 


(This  appendix  is  extracted  from  CTI-TR-9207  report  currently  in  preparation.] 


Bayesian  Data  Assimilation 

Given  a  set  of  measured  data  taken  at  periodic  intervals  over  some  specified  spatial  mesh  in  the 
atmosphere,  the  tast  of  a  Kalman  filter  is  to  integrate  these  data  with  a  known  or  modeled 
description  of  the  expected  dynamics  to  produce  a  'best'  estimate  of  the  current  atmospheric 
state. 

In  the  present  description  the  known  atmospheric  dynamics  are  specified  by  a  set  of  equations 
of  motion  describing  the  constraints  of  continuity,  energy  and  momentum  conservation.  For  the 
planetary  boundary  layer  the  flow  is  well  modeled  as  incompressible.  The  equation  for  this 
condition  constitutes  the  continuity  constraint.  In  general  both  a  momentum  and  an  energy 
constraint  are  also  needed  to  describe  the  dynamics  of  a  3  dimensional  boundary  layer.  In  the 
analysis  presented  here  we  will  consider  only  a  neutrally  stable  atmosphere  and  will  not  include 
effects  of  latent  heat  release  and  other  effects  of  buoyancy.  In  this  limit  only  the  momentum  and 
continuity  equations  are  needed  to  describe  the  flow.  Expressing  the  velocity  as  the  curt  of  a 
stream  function  automatically  satisfies  the  continuity  equation.  The  distribution  and  evolution  of 
the  source  of  this  stream  function  (the  vorticity)  in  space  represents  the  fundamental  description 
of  the  atmospheric  state. 

In  general,  each  component  of  the  vorticity  vector  forms  the  source  function  for  a  Poisson 
equation  for  the  corresponding  component  of  the  velocity  stream  function 


V2vF  =  Q 


(1) 


The  fluid  velocity  is  the  curl  of  the  solution  of  this  equation 


U  =  V  x  Y  (2) 

By  its  definition  this  velocity  field  is  incompressible.  Given  this  prescription  for  calculating  the 
velocity  from  the  vorticity,  the  momentum  equation  for  the  vorticity  describes  the  temporal 
evolution  of  the  hydrodynamic  state: 


^  =  -u*vn-n*vu+vv2n  0) 

dt 

The  first  two  terms  on  the  right  hand  side  of  Eq.(3)  characterize  the  inherent  dynamics  of  the 
modeled  inviscid  atmosphere.  They  represent  respectively  advective  transport  of  vorticity  and 
stretching  of  vorticity.  The  last  term  represents  viscous  diffusion.  In  most  of  the  present 
analyses  we  will  set  this  term  to  zero. 

In  the  Kalman  filter  approach  to  data  assimilation  we  need  two  models  :  one  for  the  dynamics  of 
the  system  being  analysed  and  one  characterizing  the  relationship  between  this  dynamical  model 
and  the  measured  data.  The  Kalman  filter  merges  these  two  models  and  generates  a  single 
equation  describing  the  evolution  of  a  best  estimate  of  the  system  state. 
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We  have  chosen  vorticity  as  the  fundamental  state  variable.  The  evolution  of  the  vorticity 
distribution  can  be  represented  in  either  of  two  modes:  Eulerian  or  Lagrangian.  In  the  Eulerian 
mode  the  state  variable  is  the  3D  vector  vorticity  evaluated  at  the  nodes  of  a  prespecified  3D 
spatial  grid.  In  the  Lagrangian  mode  discrete  vortex  elements  are  tracked  as  they  move  through 
space  with  the  flow.  Here  the  state  variable  has  six  component  values  for  each  vortex  element: 
3  vector  components  of  circulation  plus  3  position  components. 

In  the  next  section  the  Kalman  filter  formalism  will  be  developed  for  the  Eulerian  mode  and  in 
the  following  section  for  the  Lagrangian  mode.  In  both  cases  the  filter  has  the  basic  form 

^  =  F(X)  +  \  P(X,  X' )  dH(^p-{Z(S)  -  H(S, {X})}dS dX'  (4) 

Here  X  is  the  estimate  of  the  state  variable,  F(X)  the  (known  or  expected)  state  dynamics,  Z(S) 
the  measured  data  at  the  data  space  point  S,  H(S,X)  the  prediction  of  the  data  at  S  for  the  state 
variable  set  {X},  and  P(X,X')  is  the  expected  covariance  of  the  estimate  X.  P(X,X’)  is  to  be 
calculated  from  an  ancillary  equation  (the  'correct'  Kalman  filter)  or  is  estimated  empirically(sub- 
optimal  filtering).  In  the  form  of  Eq.(4)  we  have  assumed  uncorrelated,  unit  variance,  white 
Gaussian  noise  for  the  measurement.  In  this  form  Eq.  (4)  may  be  thought  of  as  being  the  result 
of  a  least  squares  fit  of  the  data  model  to  the  measurements. 


EULERIAN  MODE  -  Tracking  vorticity 

For  uncorrelated  measurement  noise  having  unit  variance,  the  basic  Kalman  filter  for  the 
evolution  of  the  estimate  of  the  vorticity  state  of  an  inviscid  fluid  has  the  general  form 


p)  O  A  A  -  M 

^  =  -U  •  V  Q-  Q.  VU  +  J  dx"  P(x,  x"  )£  J  dx‘H(x,x " ,  8) 


Z(x  ,6)-  Zprtd (x  ,  8;  Q) 


(5) 

Here  Zpnd  is  the  current  predicted  value  of  the  data  measurement  vector  based  on  the  current 

A 

estimate  of  the  vorticity  density  field  £l(x). 


Zpnd(x  ; Q,  8, )  =  | H(x  , x'" ,  8, ) Q(x'" )dx"'  (6) 

where  qj  is  the  inclination  of  the  ith  line  of  sight; 

The  last  term  in  Eq.(5)  is  the  update  correction  to  the  predicted  fluid  dynamics  that  the  Kalman 
filter  prescribes  to  account  for  deviations  between  the  actual  observations  and  the  current 
predictions. 

For  the  purpose  of  demonstrating  the  use  of  the  Kalman  filter  for  the  assimilation  of  lidar 
measurements  with  hydrodynamic  models,  we  have  limited  the  calculations  presented  in  this 
note  to  a  two  dimensional  atmosphere  having  only  height  and  a  single  horizontal  coordinate. 
However,  for  an  actual  implementation  with  real  atmospheric  data  a  full  3D  model  will  be 
needed.  Thus  three  dimensions  will  be  retained  in  the  formulation  but  the  numerical  examples 
are  two  dimensional.  The  extension  to  three  dimensions  is  straightforward  numerically  although 
the  implementation  requires  substantially  greater  computer  resources. 
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Generally  speaking,  two  dimensional  turbulent  flows  behave  qualitatively  differently  from  three 
dimensional  flows.  The  stretching  term  in  the  momentum  equation  (the  second  term  on  the  right 
hand  side  of  Eq.  (5))  allows  vorticity  to  be  compressed  and  is  responsible  for  the  characteristic 
one-sided  cascade  of  turbulent  energy  from  large  to  small  scales  that  typifies  fully  developed 
turbulence.  In  two  dimensions,  density  of  vorticity  is  conserved  and  the  long  term  turbulent 
transport  cascade  processes  differ.  However,  many  of  the  basic  short  term  structural 
characteristics  of  the  turbulence  process  are  qualitatively  similar  to  those  in  3D  flows. 


The  Measurement  Transfer  Function  Hfx'.x",  qj) 

The  measurement  function  relates  the  measured  data  to  the  function  describing  the  atmospheric 
state.  Since  we  have  selected  vorticity  as  the  fundamental  state  variable  and  radial  velocity 
along  selected  lines  of  sight  for  the  measurement,  the  transfer  function  H  includes  the  process  of 
solving  the  Poisson  equation  to  get  the  stream  function,  taking  its  curl  to  get  velocity,  forming 
the  projection  dot  product  to  get  the  radial  component,  masking  the  2D  continuous  distribution  to 
limit  the  prediction  to  include  just  the  sampled  lines  of  sight,  and  finally,  applying  the  range 
response  function  of  the  lidar.  The  first  three  of  these  steps  as  well  as  the  last  are  homogeneous 
operations  and  can  be  represented  as  spatial  convolution  operations.  The  fourth  is  equivalent  to 
multiplying  by  a  spatial  mask.  Thus  it  is  convenient  to  write  H(x',x")  in  the  form  of  a  product  of 
two  single  variable  functions;  a  homogeneous  term  H0  dependent  only  on  the  vector 
displacement  between  the  two  space  points  and  a  space  local  term  H-j  - 

H(x,x")  =  H0(x-x)H}(x')  (7) 

The  first  term  is  derived  from  the  definition  of  the  radial  velocity 


ZH(x;/)  =  n(/).Vx>P(x)  (8) 


A 

where  n(/)  is  the  unit  vector  defining  the  Ith  lidar  line  of  sight. 

As  will  become  evident,  we  will  invoke  a  spatially  homogeneous  constraint  to  limit  the  complexity 
of  the  analysis.  For  such  systems,  many  of  the  integral  relations  are  most  easily  evaluated  in  the 
spatial  Fourier  domain.  In  anticipation  of  this  we  provide  here  various  Fourier  domain  forms  of 
the  above  expressions. 

The  Fourier  transform  of  Eq.  (8) 


Z,„,(k ;/)  =  H_(k;/)  fi(k)  =  -/{n(/)  k/|k|!}xfi(k)  (9) 

Interchanging  the  dot  and  cross  operation  allows  us  to  identify  the  measurement  transfer  function 
as 

Hs(k;/)  =  -/n(/)xk/|k|2  (10) 

For  narrow  beam  lidar  sampling,  the  second  term  in  Eq.  (7)  is  a  delta  function  defining  the 
locations  sampled  by  the  lidar  beam 
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Hx(x\l)  =  Hx(x,y,zJ)  =  5(z-xcos(^)tan(0))5(2->'sin(9)  tan(0)) 
where  q  and  j  are  the  polar  and  aximuthal  angles  of  the  Ith  line  of  sight  of  the  iidar. 


(11) 


The  State  Covariance  term  P(x\x") 

The  state  covariance  term  in  a  Kalman  filter  is  supposed  to  be  evaluated  from  an  ancillary 
differential  equation  that  continually  estimates  and  updates  the  precision  and  expected 
correlations  of  the  predicted  state  estimates.  In  a  true  Kalman  implementation  this  a  non-linear 
Ricatti  type  equation  for  a  two  variable  function,  and  can  dominate  the  computational  load  when 
the  number  of  state  variable  components  is  large.  In  most  applications  to  the  estimation  of 
multidimensional  scenes,  the  computational  requirements  for  this  evaluation  are  extreme  and,  as 
a  result,  sub  optimal  prescriptive  models  for  Pfc'.x")  are  usually  used.  This  is  the  approach  we 
will  adopt  here.  A  similar  approach  is  commonly  used  in  weather  forecasting  where  the  weather 
estimates  are  generated  using  a  procedure  similar  in  principle  to  Eq.  (5)  with  empirical 
prescriptions  used  for  the  influence  function  P(x,,x").  These  algorithms  are  often  refered  to  as 
01  methods  ('Optimal  Interpolation'  )  to  distinguish  them  from  true  Bayesian  or  Kalman 
prescriptions  in  which  a  formally  derived  statistical  algorithm  is  used  to  estimate  Pfr'.x").  In  this 
note  we  will  use  the  generic  term  Kalman  filter  to  refer  to  both  forms  of  data  analysis. 

The  term  P(x',x")  term  is  formally  defined  as  the  covariance  of  the  error  of  the  state  estimate 
and  can  be  expressed  as 


p(x',x")  =<  (Q'-n')(nM-nM)  > 


(12) 


where  W  is  the  true  vorticity  at  the  point  x\  H'  its  estimate, and  the  brackets  <>  mean  expected 
value.  When  the  vorticity  is  spatially  uncorrelated  and  white,  Eq.  (5)  reduces  to  a  least  square 
estimate  (  it  actually  becomes  the  maximum  likelihood  estimate:  in  Eq.  (5)  the  measurement 
noise  has  been  assumed  white  and  uncorrelated  so  that  the  ML  result  is  simply  the  least  square 
value). 

When  the  state  estimates  are  all  uncorrelated,  each  of  the  terms  in  Eq.  (5)  can  be  easily 
interpreted.  Replacing  P(x’,x")  by  P0  times  the  delta  function  and  integrating  gives 


r)Cy  A  A  - 

~  =  -U  •  V  Q-  Q.  VU  +  P0  Y  f  dx'H(x  ,  x,  0) 

ot  a  j 


Z(x  ,  0)  -  Zpnd(x  ,  0;  Q) 


(13) 


The  integral  term  in  Eq.(13)  can  be  recognized  as  proportional  to  the  derivative  of  the  error 
energy  (integrated  over  all  space)  with  respect  to  the  vorticity  estimate  at  the  space  point  x’: 


dE 

0Q 


-°-5-4zK 

an  * 


Z (x  ,  0)  -  Zpnd  (x  ,  0;  Q) 


(14) 
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Thus  the  Bayesian  estimation  process  embedded  in  the  Kalman  filter  (Eq.  13)  adjusts  the  current 
estimate  of  the  vorticity  state  at  each  time  step  in  two  ways.  First  the  estimate  is  advected  and 
concentrated  according  to  the  normal  fluid  dynamic  rules  for  vorticity  dynamics  of  inviscid  flow. 
Second  the  entire  vorticity  field  is  altered  at  each  point  in  space  according  to  one  step  of  a 
steepest  descent  prescription  whose  goal  is  to  minimize  the  current  mean  square  discrepancy 
between  the  actual  data  and  the  current  predictions. 

This  latter  independent  adjustment  of  the  vorticity  at  each  point  in  space  is  appropriate  only 
when  the  individual  vorticity  values  are  truly  independent.  A  finite  Pfr'.x")  function  allows  for 
correlation  between  neighboring  points.  The  integration  over  x*  in  Eq  (13)  blurs  out  the 
connections  between  the  data  error  terms  and  the  individual  vorticity  state  values  and  allows  the 
vorticity  correlations  predicted  dynamically  to  be  maintained  during  the  data  update  process. 

In  a  suboptimal  analysis  we  need  to  select  an  appropriate  form  for  P(x’,xM).  For  a 
homogeneous  boundary  layer  we  can  require  that  P  be  only  a  function  of  the  (vector)  horizontal 
difference  distance  x'-x"  together  with  a  parametric  dependence  on  the  height.  To  simplify  the 
2D  demonstration  analysis  we  will  assume  initially  a  dependence  on  the  3D  vector  difference 
displacement.  This  will  allow  Fourier  transform  methods  to  be  used  to  simplify  the  computation 
of  the  update  of  the  vorticity  state.  In  this  limit  P  is  described  completely  by  its  Fourier  transform 
P(k).  In  the  present  analysis  we  will  assume  simple,  parametric  forms  for  P(k).  In  a  more 
extended  analysis  it  would  be  appropriate  to  develop  more  direct  models  that  more  accurately 
represent  the  full  Kalman  formalism. 


LAGRANGIAN  MODE  -  Tracking  Vortices 

Instead  of  treating  the  vorticity  density  as  the  fundamental  state  variable,  we  imagine  a  model  of 
the  fluid  in  which  a  cloud  of  preexistent  discrete  vortex  elements  are  the  source  of  the  velocity 
field.  In  expectation  the  density  of  this  cloud  is  the  vorticity.  However,  the  computation  tracks 
individual  vortex  elements  (their  location,  orientation,  and  strength).  The  velocity  of  any  element 
can  be  calculated  by  directly  summing  the  Green's  function  contribution  from  all  other  elements. 
Alternately  it  may  be  evaluated  by  adding  up  the  contributions  of  all  vortex  elements  in  each 
mesh  cell  to  get  a  vorticity  density,  solving  a  Poisson  equation  for  the  stream  function  throughout 
the  mesh,  and  interpolating  the  value  of  the  stream  function  curl  at  the  vortex  location  to  get  the 
velocity.  (Using  fast  Poisson  solvers,  the  first  method  requires  of  order  N2  operations  per  time 
step  where  N  is  the  number  of  vortex  elements,  whereas  the  second  requires  of  order  N+M 
log(M)  with  M  the  number  of  mesh  cells.) 

Both  Lagrangian  and  Eulerian  methods  are  commonly  used  in  computing  three  dimensional  fluid 
flows  and  both  have  their  advantages  and  disadvantages.  For  simply  computing  the  flow  of  a 
fluid  both  methods  remain  relatively  closely  related  in  their  computational  implementations. 
However,  in  the  Kalman  filter  assimilation  of  concurrent  data,  distinct  differences  arise  that  have 
conceptually  different  mathemetical  interpretations.  These  arise  primarily  from  the  manner  in 
which  the  system  state  is  defined.  In  the  Eulerian  approach,  the  three  component  vector  vorticity 
is  the  natural  state  variable,  and  the  filter  has  the  formal  form  of  an  equation  for  3D  vector  fluid 
evolving  in  3D  space  under  the  influence  of  a  distributed  source  function  generated  by  local 
differences  between  measured  and  predicted  data  values.  In  the  Lagrangian  approach,  the  state 
variable  is  a  six  component  vector  (vortex  element  circulation  and  spatial  location  (both  3D 
vectors))  and  the  formal  form  of  the  Kalman  filter  is  now  a  cloud  of  6D  vector  elements  evolving 
in  a  6D  space.  This  cloud  has  fluid-like  qualities  and  the  filter  analysis  can  be  recast  formally 
into  a  fluid  description  in  6D  state  space.  This  space  can  have  conceptually  different  properties 
than  the  3D  spatial  space  used  in  the  Eulerian  description. 
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The  Lagrangian  state  variable,  X=(a,x),  has  six  vector  components  for  each  vortex  element;  x  is 
the  position  of  the  element  and  a  its  circulation.  The  flow  field  associated  with  the  elementary 
vortex  element  located  at  xp  and  having  circulation  ap  is  described  by  its  stream  function: 

vJ,  =  ap/(x-xp)  (15) 


where  f  is  the  Green's  function  for  a  unit  source  at  the  origin: 

/(x)  =  l/|x|  (16) 

The  fluid  velocity  u  and  its  projection  H  on  the  the  lidar  line  of  sight  are  given  by 


P 

(17) 

//(x.x^a,)  =  n •  u  =  £n •  ap  x  V/(x - xp)  =  £ n  x  ap  •  V/(x- xp) 

p  p 

The  Lagrangian  update  equation  for  the  Kalman  filter  is  derived  from  Eq.(4).  To  simplify  the 
analysis  we  will  assume  a  particular  form  for  the  covariance  P(X,X’),  i.e., 


P(X,X,)=P((x,a).(x\a'))= 


'P(x-x') 
L  0 


0  ' 
o^5(a-a')> 


(18) 


The  first  assumption  that  P  depends  only  on  the  separation  in  the  6D  state  space  between  the 
two  points  (x.a)  and  (x'.a')  is  equivalent  to  assuming  homogeneity  (no  preferred  spatial  origin,  no 
preferred  vortex  orientation  or  strength  and  should  be  valid  in  homgeneous  turbulence.  In  a 
boundary  layer,  however,  the  assumption  fails  for  turbulence  scales  comparable  to  the  local 
distance  from  the  bounding  surface). 

The  second  assumption  that  the  circulation  estimates  are  delta  function  correlated  is  introduced 
without  justification  and  is  purely  for  convenience  ~f  *he  initial  analysis. 

In  this  approximation  the  update  equation  for  the  vortex  element  positions  becomes  (from  Eq. 
(4)) 


^  =  F<*,,a,)+ j  />(*,  -  ,.)ag(‘^,,,){zW  *  H(s,{X))}dsdx'  (19) 
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Because  of  the  homogeneity  assumption  we  can  rewrite  the  data  term  in  Eq(19)  as  the  gradient 
of  a  psuedo  potential  evaluated  at  the  location  of  the  vortex: 


d\p 

dt 


F(x,,.,)+vr(*,.) |„w 


(20) 


where 

F(x,a)  =  |  P(x-\,)H(six\^){Z(s)-H(s,{X})}dsdx'  (21) 


Eqs.  (20)  and  (21)  demonstrate  very  clearly  the  basic  mechanics  of  the  Kalman  filter  when 
applied  to  distributed  measurements  that  sparsely  sample  a  continuous  dynamical  field.  As  a 
result  of  the  basic  Bayesian  formulation  we  expect  (i.e.,  we  demanded)  that  the  dynamical  model 
and  the  data  be  combined  to  yield  a  'best'  estimate  of  the  current  state.  This  best  estimate 
criterion  is  defined  as  being  that  which  produces  the  minimum  expected  uncertainty  of  the 
estimated  state  (minimum  variance  in  the  present  description  because  of  the  assumption  of 
uncorrelated  measurement  noise  in  Eq.  (4)).  When  the  measured  data  disagree  with  current 
predictions  Eq.  (21)  establishes  an  error  signal  at  those  locations  of  the  disagreement.  This  error 
signal  is  multiplied  by  the  expected  value  of  the  measurement  to  form  an  error  source  function. 
This  latter  multiplication  is  necessary  so  that  finite  measurement  values  originating  from  regions 
from  where  no  measured  levels  are  expected  are  properly  interpreted  as  being  noise  and  will  be 
rejected.  This  source  signal  is  then  convolved  with  the  state  covariance  model  P  to  create  a 
scalar  potential  field  V.  This  latter  convolution  ensures  that  all  state  points  in  the  vicinity  of  an 
error  signal  that  are  correlated  will  respond  to  the  error.  The  gradient  of  this  potential  field 
estabishes  an  advection  velocity  field  that  will  induce  a  motion  designed  to  reduce  the  error 
signal.  This  velocity  is  added  to  the  physical  velocity  field  due  to  the  interaction  between  all  the 
currently  estimated  vortices. 

The  advantage  of  the  convolution  form  for  establishing  the  error  advection  field  cannot  be 
understated.  A  major  problem  in  steepest  descent  relaxation  methods  (of  which  this  in  one)  is 
the  difficulty  of  getting  the  whole  solution  field  to  respond  to  highly  localized  errors  without 
getting  stuck  in  local  minima.  By  having  the  error  signals  blurred  out  by  the  covariance  function, 
it  is  possible  to  make  all  points  in  the  field  feel  some  influence  of  every  error  signal,  local  as  well 
as  distant,  and  thus  allow  all  points  to  respond  instantaneously  and  continuously.  Because  the 
response  at  a  given  vortex  depends  on  the  coordinates  and  strengths  of  all  other  vortices,  a 
smooth  type  of  response  can  be  expected  if  the  effective  psuedo  force  field  is  suitably  chosen. 

In  principle,  the  influence  function  P  should  be  calculated  from  a  separate  covariance  evolution 
equation.  Without  the  homogeneity  assumption  to  limit  the  functional  dependence  this  can  be  a 
horrendously  difficult  task.  For  homogeneous  environments,  P  is  simply  a  scalar  function  of  one 
space  variable  and  should  be  easily  computed  (especially  since  it  depends  only  on  the  expected 
data  statistics,  not  the  actual  data  itself).  In  the  present  analysis,  however,  we  do  not  attempt  to 
solve  directly  for  P  but  will  arbitrarily  prescribe  forms  heuristically. 

Substituting  for  H  from  Eq.(  17)  and  setting  F(x)  equal  to  the  local  fluid  velocity  gives 

^  =  u(x)  +  v[J  P(x  -  x')tf(s,*'){Z(s)-  ff(s,{X})}<*<*[’]  (22) 


In  Eq.  (22)  the  integration  over  s  is  limited  to  the  actual  measured  points.  For  any  measurement  system 
which  simply  samples  points  in  the  state  space  it  is  convenient  to  introduce  a  mask  function  Hs(x,n)  to 
accomplish  this  constraint.  Invoking  also  the  homogeneous  condition  for  H  then  gives 
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~  =  u(x)  + V^(x,a) 
at 


(23) 


where  V  is  the  scalar  'pseudo-potential'  defined  by 


V(x,„  )  =  jp(i-x')£ 


H (x,  -  X*  )H,  (x„  n)\  Z(x, )  -  £  tf(x,  -  xq )  Wx,  k*  (24) 


where  now  the  integrations  over  x*  and  xs  extend  over  the  full  state  space.  In  recognition  of  the 
fact  that  H  depends  on  the  lidar  look  direction  n,  we  have  separated  the  contributions  from 
separate  look  directions  into  a  sum  over  these  directions. 

Equations  23  and  24  are  convenient  for  numerical  integrations  using  Fourier  methods  for 
evaluating  the  advection  fields.  Using  the  definitions  for  u  and  H 

u(*,)  =  XapxV/(x-xp)  (25) 

p 

H(5x, a, n)  =  n •  a  x  Vf  ( 5x )  =  nxa-V/ (5x) 


we  may  write  the  psuedo-potential  function  in  the  form 

V(x, a)  =  J  P(x  -  x' )£ [n  x  a  •  V/ (x,  - x’ )Z,(x, , n)]c/x’ ds  (26) 

n 

where  Zs(x,  n)  is  the  (masked)  error  source  function  for  the  line  of  sight  n: 


ZJ(x,n)  =  //J(x,n)|Z(x)-^//(x-x?,n,aq) 

The  Fourier  transforms  of  H  and  of  the  true  velocity  fields  are 

H(  k,  n,  a)  =  /'n  x  a  •  l/(k) 

u(k)  =  /'Z  aP  ><H/'(k)exp(/k-xp) 

p 

The  psuedo-potential  may  be  written  in  the  form 

F(k,a)  =  -/P(k)/(k)k-a  xZ,(k)/7?AsAxA/ 


(27) 


(28) 


(29) 


where  we  have  defined  the  masked  measurement  error  vector  Zs(x)  as  a  vector  sum  over  all 
lines  of  sight 
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Z,(x)  =  £nZ,(x,n) 


(30) 


This  quantity  is  simply  the  3D  distribution  of  sensed  error  radial  velocity  (converted  to  radial 
vectors  at  each  measurement  point  and  set  to  zero  at  all  other  points). 

Equation(29)  can  be  rewritten  in  terms  of  the  curl  of  this  quantity,  i.e.,  in  terms  of  a  measurement 
error  vorticity  field  0,(x)  =  V  xZ,(x): 

V(k,  a)  =  P(  k) /(k)a  •  Cl,  (k)  /  RAsAxAt  (31) 

Thus  the  scalar  pseudo-potential  V  is  equal  to  the  projection  of  the  vortex  circulation  vector  on  to 
a  vector  potential  function  quantity  A,(x): 

F(x,a)  =  a-A,(x)  (32) 

This  vector  potential  function  can  be  interpreted  as  being  a  vector  field  generated  by  the 
(vector)  source  function  ^(x)  with  the  covariance  function  P(x)  serving  as  the  Green's 
function.  In  Fourier  space 

A,(k)  =  PCk^Ck)  /  RAxAsAt.  (33) 


In  turn  the  qi-antity  ^(x)  is  the  stream  function  associated  with  the  pseudo  vorticity  error  field 
fls  which  is  derived  from  the  observed  line-of-  sight  velocities  according  to  Eqs.  27  and  30: 


y,(k)=/(k)n,(k) 


(34) 


In  other  words 


4(k)  =  P(k)/(k)n,(k)/i?AsAxA/  (35) 


Equations  (23)  to  (34)  provide  a  means  for  calculating  the  evolution  of  positions  of  the  discrete 
vortex  elements  given  the  input  data  and  the  values  of  the  vortex  strengths. 

Equations  describing  the  evolution  of  the  circulations  of  the  vortex  elements  may  be  constructed 
using  a  similar  analysis.  From  Eq.  (4)  and  (18),  and  using  the  same  assumptions  used  in 
deriving  Eq.  (22)  to  (30),  we  obtain  an  expression  for  the  evolution  of  the  estimates  of  the  vortex 
element  strengths: 


dt 


=  J  gffia  - a„)  dZp(Xp  ^,,a’n^Zs(x,)dadx,  /  flAsAf 


(36) 
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Eq  (25)  allows  this  equation  to  be  rewritten  in  the  form 


^  =  Xn  x  j crS(a  - ap)V/(xp -x,)Z,(x„n)dx,da/ RAsAt  (37) 

Ul  n 

Expressed  in  the  Fourier  domain 

— -  =  -i cra  [ kxZ,(k)/(k)exp(/k*x  )dk/ RAsAt 
at  J 

=  -c^j /(k)Q,(k)exp(/k-xp)tfk/i?AsA/  (38) 

In  other  words,  the  growth  rate  of  the  strength  of  vortex  element  p  is  just  proportional  to  the 
negative  of  the  pseudo-stream  function  evaluated  at  the  vortex  location 


d*  o^FfxJ 

p  _  a  jv  p/ 

~dt~  RASAt 


(39) 
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